text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	storage(EP::Model, inputs::Dict, Reserves::Int, OperationWrapping::Int, LongDurationStorage::Int, EnergyShareRequirement::Int, CapacityReserveMargin::Int)
A wide range of energy storage devices (all $o \in \mathcal{O}$) can be modeled in GenX, using one of two generic storage formulations: (1) storage technologies with symmetric charge and discharge capacity (all $o \in \mathcal{O}^{sym}$), such as Lithium-ion batteries and most other electrochemical storage devices that use the same components for both charge and discharge; and (2) storage technologies that employ distinct and potentially asymmetric charge and discharge capacities (all $o \in \mathcal{O}^{asym}$), such as most thermal storage technologies or hydrogen electrolysis/storage/fuel cell or combustion turbine systems.
**Storage with symmetric charge and discharge capacity**
For storage technologies with symmetric charge and discharge capacity (all $o \in \mathcal{O}^{sym}$), charge rate, $\Pi_{o,z,t}$, is constrained by the total installed power capacity, $\Omega_{o,z}$. Since storage resources generally represent a `cluster' of multiple similar storage devices of the same type/cost in the same zone, GenX permits storage resources to simultaneously charge and discharge (as some units could be charging while others discharge), with the simultaenous sum of charge, $\Pi_{o,z,t}$, and discharge, $\Theta_{o,z,t}$, also limited by the total installed power capacity, $\Delta^{total}_{o,z}$. These two constraints are as follows:
```math
\begin{aligned}
	&  \Pi_{o,z,t} \leq \Delta^{total}_{o,z} & \quad \forall o \in \mathcal{O}^{sym}, z \in \mathcal{Z}, t \in \mathcal{T}\\
	&  \Pi_{o,z,t} + \Theta_{o,z,t} \leq \Delta^{total}_{o,z} & \quad \forall o \in \mathcal{O}^{sym}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
These constraints are created with the function ```storage_symmetric()``` in ```storage_symmetric.jl```.
If reserves are modeled, the following two constraints replace those above:
```math
\begin{aligned}
&  \Pi_{o,z,t} + f^{charge}_{o,z,t} \leq \Delta^{total}_{o,z} & \quad \forall o \in \mathcal{O}^{sym}, z \in \mathcal{Z}, t \in \mathcal{T}\\
&  \Pi_{o,z,t} + f^{charge}_{o,z,t} + \Theta_{o,z,t} + f^{discharge}_{o,z,t} + r^{discharge}_{o,z,t} \leq \Delta^{total}_{o,z} & \quad \forall o \in \mathcal{O}^{sym}, z \in \mathcal{Z}, t \in \mathcal{T} \\
\end{aligned}
```
where $f^{charge}_{o,z,t}$ is the contribution of storage resources to frequency regulation while charging, $f^{discharge}_{o,z,t}$ is the contribution of storage resources to frequency regulation while discharging, and $r^{discharge}_{o,z,t}$ is the contribution of storage resources to upward reserves while discharging. Note that as storage resources can contribute to regulation and reserves while either charging or discharging, the proxy variables $f^{charge}_{o,z,t}, f^{discharge}_{o,z,t}$ and $r^{charge}_{o,z,t}, r^{discharge}_{o,z,t}$ are created for storage resources where the total contribution to regulation and reserves, $f_{o,z,t}, r_{o,z,t}$ is the sum of the proxy variables.
These constraints are created with the function ```storage_symmetric_reserves()``` in ```storage_symmetric.jl```.
**Storage with asymmetric charge and discharge capacity**
For storage technologies with asymmetric charge and discharge capacities (all $o \in \mathcal{O}^{asym}$), charge rate, $\Pi_{o,z,t}$, is constrained by the total installed charge capacity, $\Delta^{total, charge}_{o,z}$, as follows:
```math
\begin{aligned}
	&  \Pi_{o,z,t} \leq \Delta^{total, charge}_{o,z} & \quad \forall o \in \mathcal{O}^{asym}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
These constraints are created with the function ```storage_asymmetric()``` in ```storage_asymmetric.jl```.
If reserves are modeled, the above constraint is replaced by the following:
```math
\begin{aligned}
	&  \Pi_{o,z,t} + f^{charge}_{o,z,t} \leq \Delta^{total, charge}_{o,z} & \quad \forall o \in \mathcal{O}^{asym}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
where $f^{+}_{y=o,z,t}$ is the contribution of storage resources to frequency regulation while charging.
These constraints are created with the function ```storage_asymmetric_reserves()``` in ```storage_asymmetric.jl```.
**All storage resources**
The following constraints apply to all storage resources, $o \in \mathcal{O}$, regardless of whether the charge/discharge capacities are symmetric or asymmetric.
The following two constraints track the state of charge of the storage resources at the end of each time period, relating the volume of energy stored at the end of the time period, $\Gamma_{o,z,t}$, to the state of charge at the end of the prior time period, $\Gamma_{o,z,t-1}$, the charge and discharge decisions in the current time period, $\Pi_{o,z,t}, \Theta_{o,z,t}$, and the self discharge rate for the storage resource (if any), $\eta_{o,z}^{loss}$.  The first of these two constraints enforces storage inventory balance for interior time steps $(t \in \mathcal{T}^{interior})$, while the second enforces storage balance constraint for the initial time step $(t \in \mathcal{T}^{start})$.
```math
\begin{aligned}
	&  \Gamma_{o,z,t} =\Gamma_{o,z,t-1} - \frac{1}{\eta_{o,z}^{discharge}}\Theta_{o,z,t} + \eta_{o,z}^{charge}\Pi_{o,z,t} - \eta_{o,z}^{loss}\Gamma_{o,z,t-1}  \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}^{interior}\\
	&  \Gamma_{o,z,t} =\Gamma_{o,z,t+\tau^{period}-1} - \frac{1}{\eta_{o,z}^{discharge}}\Theta_{o,z,t} + \eta_{o,z}^{charge}\Pi_{o,z,t} - \eta_{o,z}^{loss}\Gamma_{o,z,t+\tau^{period}-1}  \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}^{start}
\end{aligned}
```
When modeling the entire year as a single chronological period with total number of time steps of $\tau^{period}$, storage inventory in the first time step is linked to storage inventory at the last time step of the period representing the year. Alternatively, when modeling the entire year with multiple representative periods, this constraint relates storage inventory in the first timestep of the representative period with the inventory at the last time step of the representative period, where each representative period is made of $\tau^{period}$ time steps. In this implementation, energy exchange between representative periods is not permitted. When modeling representative time periods, GenX enables modeling of long duration energy storage which tracks state of charge between representative periods enable energy to be moved throughout the year. If ```LongDurationStorage=1``` and ```OperationWrapping=1```, this function calls ```long_duration_storage()``` in ```long_duration_storage.jl``` to enable this feature.
The next constraint limits the volume of energy stored at any time, $\Gamma_{o,z,t}$, to be less than the installed energy storage capacity, $\Delta^{total, energy}_{o,z}$. Finally, the maximum discharge rate for storage resources, $\Pi_{o,z,t}$, is constrained to be less than the discharge power capacity, $\Omega_{o,z,t}$ or the state of charge at the end of the last period, $\Gamma_{o,z,t-1}$, whichever is lessor.
```math
\begin{aligned}
	&  \Gamma_{o,z,t} \leq \Delta^{total, energy}_{o,z} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}\\
	&  \Theta_{o,z,t} \leq \Delta^{total}_{o,z} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}\\
	&  \Theta_{o,z,t} \leq \Gamma_{o,z,t-1} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
The above constraints are established in ```storage_all()``` in ```storage_all.jl```.
If reserves are modeled, two pairs of proxy variables $f^{charge}_{o,z,t}, f^{discharge}_{o,z,t}$ and $r^{charge}_{o,z,t}, r^{discharge}_{o,z,t}$ are created for storage resources, to denote the contribution of storage resources to regulation or reserves while charging or discharging, respectively. The total contribution to regulation and reserves, $f_{o,z,t}, r_{o,z,t}$ is then the sum of the proxy variables:
```math
\begin{aligned}
	&  f_{o,z,t} = f^{charge}_{o,z,t} + f^{dicharge}_{o,z,t} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}\\
	&  r_{o,z,t} = r^{charge}_{o,z,t} + r^{dicharge}_{o,z,t} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
The total storage contribution to frequency regulation ($f_{o,z,t}$) and reserves ($r_{o,z,t}$) are each limited specified fraction of installed discharge power capacity ($\upsilon^{reg}_{y,z}, \upsilon^{rsv}_{y,z}$), reflecting the maximum ramp rate for the storage resource in whatever time interval defines the requisite response time for the regulation or reserve products (e.g., 5 mins or 15 mins or 30 mins). These response times differ by system operator and reserve product, and so the user should define these parameters in a self-consistent way for whatever system context they are modeling.
```math
\begin{aligned}
	f_{y,z,t} \leq \upsilon^{reg}_{y,z} \times \Delta^{total}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{W}, z \in \mathcal{Z}, t \in \mathcal{T} \\
	r_{y,z, t} \leq \upsilon^{rsv}_{y,z}\times \Delta^{total}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{W}, z \in \mathcal{Z}, t \in \mathcal{T}
	\end{aligned}
```
When charging, reducing the charge rate is contributing to upwards reserve and frequency regulation as it drops net demand. As such, the sum of the charge rate plus contribution to regulation and reserves up must be greater than zero. Additionally, the discharge rate plus the contribution to regulation must be greater than zero.
```math
\begin{aligned}
	&  \Pi_{o,z,t} - f^{charge}_{o,z,t} - r^{charge}_{o,z,t} \geq 0 & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}\\
	&  \Theta_{o,z,t} - f^{discharge}_{o,z,t} \geq 0 & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
Additionally, when reserves are modeled, the maximum charge rate and contribution to regulation while charging can be no greater than the available energy storage capacity, or the difference between the total energy storage capacity, $\Delta^{total, energy}_{o,z}$, and the state of charge at the end of the previous time period, $\Gamma_{o,z,t-1}$. Note that for storage to contribute to reserves down while charging, the storage device must be capable of increasing the charge rate (which increase net load).
```math
\begin{aligned}
	&  \Pi_{o,z,t} + f^{charge}_{o,z,t} \leq \Delta^{energy, total}_{o,z} - \Gamma_{o,z,t-1} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
Finally, the constraints on maximum discharge rate are replaced by the following, to account for capacity contributed to regulation and reserves:
```math
\begin{aligned}
	&  \Theta_{o,z,t} + f^{discharge}_{o,z,t} + r^{discharge}_{o,z,t} \leq \Delta^{total}_{o,z} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}\\
	&  \Theta_{o,z,t} + f^{discharge}_{o,z,t} + r^{discharge}_{o,z,t} \leq \Gamma_{o,z,t-1} & \quad \forall o \in \mathcal{O}, z \in \mathcal{Z}, t \in \mathcal{T}
\end{aligned}
```
The above reserve related constraints are established by ```storage_all_reserves()``` in ```storage_all.jl```
"""
function storage(EP::Model, inputs::Dict, Reserves::Int, OperationWrapping::Int, LongDurationStorage::Int, EnergyShareRequirement::Int, CapacityReserveMargin::Int, StorageLosses::Int)

	println("Storage Resources Module")

	dfGen = inputs["dfGen"]
	G = inputs["G"]
	T = inputs["T"]
	STOR_ALL = inputs["STOR_ALL"]

	if !isempty(STOR_ALL)
		EP = investment_energy(EP, inputs)
		EP = storage_all(EP, inputs, Reserves, OperationWrapping)

		# Include LongDurationStorage only when modeling representative periods and long-duration storage
		if OperationWrapping == 1 && !isempty(inputs["STOR_LONG_DURATION"])
			EP = long_duration_storage(EP, inputs)
		end
	end

	if !isempty(inputs["STOR_ASYMMETRIC"])
		EP = investment_charge(EP, inputs)
		EP = storage_asymmetric(EP, inputs, Reserves)
	end

	if !isempty(inputs["STOR_SYMMETRIC"])
		EP = storage_symmetric(EP, inputs, Reserves)
	end

	# ESR Lossses
	if EnergyShareRequirement >= 1
		@expression(EP, eESRStor[ESR=1:inputs["nESR"]], sum(inputs["dfESR"][:,ESR][z]*StorageLosses*sum(EP[:eELOSS][y] for y in intersect(dfGen[dfGen.Zone.==z,:R_ID],STOR_ALL)) for z=findall(x->x>0,inputs["dfESR"][:,ESR])))
		EP[:eESR] -= eESRStor
	end

	# Capacity Reserves Margin policy
	if CapacityReserveMargin > 0
		@expression(EP, eCapResMarBalanceStor[res=1:inputs["NCapacityReserveMargin"], t=1:T], sum(dfGen[y,Symbol("CapRes_$res")] * (EP[:vP][y,t] - EP[:vCHARGE][y,t])  for y in STOR_ALL))
		EP[:eCapResMarBalance] += eCapResMarBalanceStor
	end

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
"""Read REISE input matfiles, return parsed relevant data in a Dict."""
function read_case(filepath)
    println("Reading from folder: " * filepath)

    println("...loading case.mat")
    case_mat_file = MAT.matopen(joinpath(filepath, "case.mat"))
    mpc = read(case_mat_file, "mpc")

    # New case.mat analog
    case = Dict()

    # AC branches
    # dropdims() will remove extraneous dimension
    case["branchid"] = dropdims(mpc["branchid"]; dims=2)
    # convert() will convert float array to int array
    case["branch_from"] = convert(Array{Int,1}, mpc["branch"][:, 1])
    case["branch_to"] = convert(Array{Int,1}, mpc["branch"][:, 2])
    case["branch_reactance"] = mpc["branch"][:, 4]
    case["branch_rating"] = mpc["branch"][:, 6]

    # DC branches
    if "dcline" in keys(mpc)
        if isa(mpc["dclineid"], Int)
            case["dclineid"] = Int64[mpc["dclineid"]]
        else
            case["dclineid"] = dropdims(mpc["dclineid"]; dims=2)
        end
        case["dcline_from"] = convert(Array{Int,1}, mpc["dcline"][:, 1])
        case["dcline_to"] = convert(Array{Int,1}, mpc["dcline"][:, 2])
        case["dcline_pmin"] = mpc["dcline"][:, 10]
        case["dcline_pmax"] = mpc["dcline"][:, 11]
    else
        case["dclineid"] = Int64[]
        case["dcline_from"] = Int64[]
        case["dcline_to"] = Int64[]
        case["dcline_pmin"] = Float64[]
        case["dcline_pmax"] = Float64[]
    end

    # Buses
    case["busid"] = convert(Array{Int,1}, mpc["bus"][:, 1])
    case["bus_demand"] = mpc["bus"][:, 3]
    case["bus_zone"] = convert(Array{Int,1}, mpc["bus"][:, 7])

    # Generators
    case["genid"] = dropdims(mpc["genid"]; dims=2)
    genfuel = dropdims(mpc["genfuel"]; dims=2)
    case["genfuel"] = convert(Array{String,1}, genfuel)
    case["gen_bus"] = convert(Array{Int,1}, mpc["gen"][:, 1])
    case["gen_status"] = mpc["gen"][:, 8]
    case["gen_pmax"] = mpc["gen"][:, 9]
    case["gen_pmin"] = mpc["gen"][:, 10]
    case["gen_ramp30"] = mpc["gen"][:, 19]

    # Generator costs
    case["gencost"] = mpc["gencost"]

    # Load all relevant profile data from CSV files
    println("...loading demand.csv")
    case["demand"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "demand.csv")))

    println("...loading hydro.csv")
    case["hydro"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "hydro.csv")))

    println("...loading wind.csv")
    case["wind"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "wind.csv")))

    println("...loading solar.csv")
    case["solar"] = DataFrames.DataFrame(CSV.File(joinpath(filepath, "solar.csv")))

    return case
end

"""Read input matfile (if present), return parsed data in a Storage struct."""
function read_storage(filepath)::Storage
    # Fallback dataframe, in case there's no case_storage.mat file
    storage = Dict(
        "enabled" => false, "gen" => zeros(0, 21), "sd_table" => DataFrames.DataFrame()
    )
    try
        case_storage_file = MAT.matopen(joinpath(filepath, "case_storage.mat"))
        storage_mat_data = read(case_storage_file, "storage")
        println("...loading case_storage.mat")
        # Convert N x 1 array of strings into 1D array of Symbols (length N)
        column_symbols = Symbol.(vec(storage_mat_data["sd_table"]["colnames"]))
        storage = Dict(
            "enabled" => true,
            "gen" => storage_mat_data["gen"],
            "sd_table" =>
                DataFrames.DataFrame(storage_mat_data["sd_table"]["data"], column_symbols),
        )
    catch e
        println("File case_storage.mat not found in " * filepath)
    end

    # Convert Dict to NamedTuple
    storage = (; (Symbol(k) => v for (k, v) in storage)...)
    # Convert NamedTuple to Storage
    storage = Storage(; storage...)

    return storage
end

"""
    read_demand_flexibility(filepath, interval)

Load demand flexibility profiles and parameters from .csv files and return them in a
DemandFlexibility struct.
"""
function read_demand_flexibility(filepath, interval)::DemandFlexibility
    # Initialize demand flexibility
    demand_flexibility = Dict(
        "duration" => interval,
        "enabled" => "not_specified",
        "interval_balance" => true,
        "rolling_balance" => true,
    )

    # Try loading the demand flexibility parameters
    demand_flexibility_parameters = DataFrames.DataFrame()
    try
        demand_flexibility_parameters = DataFrames.DataFrame(
            CSV.File(joinpath(filepath, "demand_flexibility_parameters.csv"))
        )
        println("...loading demand flexibility parameters")

        # Create a dictionary to hold the warning messages relevant to loading the
        # demand flexibility parameters
        demand_flexibility_params_warns = Dict(
            "duration" => (
                "The demand flexibility duration parameter is not defined. Will " *
                "default to being the size of the interval."
            ),
            "enabled" => (
                "The parameter that indicates if demand flexibility is enabled is not " *
                "defined. Will default to being enabled."
            ),
            "interval_balance" => (
                "The parameter that indicates if the interval load balance constraint " *
                "is enabled is not defined. Will default to being enabled."
            ),
            "rolling_balance" => (
                "The parameter that indicates if the rolling load balance constraint " *
                "is enabled is not defined. Will default to being enabled."
            ),
        )

        # Try assigning the different demand flexibility parameters from the file
        for k in keys(demand_flexibility_params_warns)
            try
                demand_flexibility[k] = demand_flexibility_parameters[1, k]
            catch e
                println(demand_flexibility_params_warns[k])
            end
        end

    catch e
        println("Demand flexibility parameters not found in " * filepath)
        println(
            "Demand flexibility parameters will default to allowing demand flexibility " *
            "to occur.",
        )
    end

    # Check the feasibility of the duration parameter
    if demand_flexibility["duration"] > interval
        @warn (
            "Demand flexibility durations greater than the interval length are " *
            "set equal to the interval length."
        )
        demand_flexibility["duration"] = interval
    end

    # Prevent the rolling_balance constraint according to the duration parameter
    demand_flexibility["rolling_balance"] &= !(demand_flexibility["duration"] == interval)

    # Try loading the demand flexibility and demand flexibility cost profiles
    for s in ["up", "dn"]
        # Pre-specify the demand flexibility and demand flexibility cost profiles
        demand_flexibility["flex_amt_" * s] = nothing
        demand_flexibility["cost_" * s] = nothing

        # Only try loading the profiles if demand flexibility is enabled
        if demand_flexibility["enabled"] == "not_specified" ||
            (demand_flexibility["enabled"])
            # Try loading the demand flexibility profiles
            try
                demand_flexibility["flex_amt_" * s] = DataFrames.DataFrame(
                    CSV.File(joinpath(filepath, "demand_flexibility_" * s * ".csv"))
                )
                println("...loading demand flexibility " * s * " profiles")
            catch e
                println("Demand flexibility " * s * " profile not found in " * filepath)
            end

            # Try loading the demand flexibility cost profiles
            try
                demand_flexibility["cost_" * s] = DataFrames.DataFrame(
                    CSV.File(joinpath(filepath, "demand_flexibility_cost_" * s * ".csv"))
                )
                println("...loading demand flexibility " * s * "-shift cost profiles")
            catch e
                println(
                    "Demand flexibility " *
                    s *
                    "-shift cost profiles not found in " *
                    filepath *
                    ". Will default to no cost for " *
                    s *
                    "-shifting demand.",
                )
            end
        end
    end

    # If demand flexibility is enabled but at least one demand flexibility profile is nothing
    if demand_flexibility["enabled"] == true && (
        isnothing(demand_flexibility["flex_amt_up"]) ||
        (isnothing(demand_flexibility["flex_amt_dn"]))
    )
        throw(
            ErrorException(
                "Demand flexibility was specified to be enabled, however at " *
                "least one demand flexibility profile is missing. Please make sure both " *
                "demand flexibility profiles are included in " *
                filepath,
            ),
        )
    elseif demand_flexibility["enabled"] == "not_specified"
        if !isnothing(demand_flexibility["flex_amt_up"]) &&
            (!isnothing(demand_flexibility["flex_amt_dn"]))
            demand_flexibility["enabled"] = true
        else
            if !isnothing(demand_flexibility["flex_amt_up"]) ||
                (!isnothing(demand_flexibility["flex_amt_dn"]))
                @warn (
                    "The exclusion of one of the demand flexibility profiles has resulted " *
                    "in demand flexibility not being enabled."
                )
            end
            demand_flexibility["enabled"] = false
        end
    end

    # Set the demand flexibility constraints to false if enabled is false
    if !demand_flexibility["enabled"]
        demand_flexibility["interval_balance"] = false
        demand_flexibility["rolling_balance"] = false
    end

    # Convert Dict to NamedTuple
    demand_flexibility = (; (Symbol(k) => v for (k, v) in demand_flexibility)...)

    # Convert NamedTuple to DemandFlexibility object
    demand_flexibility = DemandFlexibility(; demand_flexibility...)

    return demand_flexibility
end

"""
    _make_bus_demand_weighting(case)

Given a Case object, build a sparse matrix that indicates the weighting of each bus in
    each zone.
"""
function _make_bus_demand_weighting(case::Case)::SparseMatrixCSC
    bus_idx = 1:length(case.busid)
    bus_df = DataFrames.DataFrame(;
        name=case.busid, load=case.bus_demand, zone=case.bus_zone
    )
    zone_demand = DataFrames.combine(DataFrames.groupby(bus_df, :zone), :load => sum)
    zone_list = sort(collect(Set(case.bus_zone)))
    zone_idx = 1:length(zone_list)
    zone_id2idx = Dict(zone_list .=> zone_idx)
    bus_df_with_zone_load = DataFrames.innerjoin(bus_df, zone_demand; on=:zone)
    bus_share = bus_df[:, :load] ./ bus_df_with_zone_load[:, :load_sum]
    bus_zone_idx = Int64[zone_id2idx[z] for z in case.bus_zone]
    zone_to_bus_shares = sparse(bus_zone_idx, bus_idx, bus_share)::SparseMatrixCSC
    return zone_to_bus_shares
end

"""
    reformat_demand_flexibility_input(case, demand_flexibility, sets)

Inspect the raw input files of flexiblility amount and cost, convert the zone/bus mixed
    data to bus data
"""
function reformat_demand_flexibility_input(
    case::Case, demand_flexibility::DemandFlexibility, sets::Sets
)

    # check consistency of flexibility input headers
    if !all(
        sort(names(demand_flexibility.flex_amt_up)) .==
        sort(names(demand_flexibility.flex_amt_dn)),
    )
        throw(
            ErrorException(
                "The flexible bus/load zone specified in the up/down " *
                "input csvs do not match. Please check the input files to make sure " *
                "every flexible bus or load zone has corresponding columns in both " *
                "flexibility csvs.",
            ),
        )
    end

    # list of zones in the network
    zone_list = sort(collect(Set(case.bus_zone)))
    # distribute zone-level aggregated number to buses
    zone_to_bus_shares = _make_bus_demand_weighting(case)
    # incidence matrix of mapping between zone and bus
    zone_to_bus_incidence = deepcopy(zone_to_bus_shares)
    (x, y, v) = findnz(zone_to_bus_incidence)
    for i in 1:nnz(zone_to_bus_incidence)
        # buses with no load still show up as nonzero entries
        if v[i] > 0
            zone_to_bus_incidence[x[i], y[i]] = 1
        end
    end

    # create an empty demand flexibility object and copy unchanged fields
    demand_flexibility_updated = Dict(
        "duration" => demand_flexibility.duration,
        "enabled" => demand_flexibility.enabled,
        "interval_balance" => demand_flexibility.interval_balance,
        "rolling_balance" => demand_flexibility.rolling_balance,
        "flex_amt_up" => nothing,
        "flex_amt_dn" => nothing,
        "cost_up" => nothing,
        "cost_dn" => nothing,
    )

    # iterate through fields
    for field in ["flex_amt_up", "flex_amt_dn", "cost_up", "cost_dn"]
        # dataframe in the corresponding field
        fh = getfield(demand_flexibility, Symbol(field))

        # skip if no cost
        if isnothing(fh)
            continue
        end

        # headers
        flexible_str = names(fh)[2:end]

        # index of columns specifing the flexibility of a zone
        zone_columns_idx = findall(x -> occursin("zone.", x), flexible_str)

        # numeric ID corresponding to zone columns
        flexible_zone_id = [
            parse(Int64, replace(flexible_str[i], "zone." => "")) for i in zone_columns_idx
        ]
        flexible_zone_num = length(flexible_zone_id)

        # do zone-bus conversion only when zone columns are present
        if flexible_zone_num > 0

            # index of columns specifing the flexibility of a bus
            bus_columns_idx = findall(x -> !occursin(".", x), flexible_str)

            # numeric ID corresponding to bus columns
            flexible_bus_id = [parse(Int64, flexible_str[i]) for i in bus_columns_idx]
            flexible_bus_idx = [sets.bus_id2idx[x] for x in flexible_bus_id]
            flexible_bus_num = length(flexible_bus_id)

            # remove bus columns from zone to bus mapping
            zone_to_bus_shares[:, flexible_bus_idx] .= 0

            # re-normalize the rows to distribute flexiblity among un-specified buses
            for i in 1:length(zone_list)
                zone_to_bus_shares[i, :] ./= sum(zone_to_bus_shares[i, :])
            end

            # check if zone numbers are correct
            if !all([issubset(i, zone_list) for i in flexible_zone_id])
                throw(
                    ErrorException(
                        "Invalid load zone numeric ID(s) in demand flexibility input files!"
                    ),
                )
            elseif !all([issubset(i, case.busid) for i in flexible_bus_id])
                throw(
                    ErrorException(
                        "Invalid load bus numeric ID(s) in demand flexibility input files!"
                    ),
                )
            end

            # list index of each zone in input file in the sorted list of zones
            zone_cols_idx = [findfirst(y -> y == x, zone_list) for x in flexible_zone_id]

            # for flex amt, the zone numbers are the total flexibility in the zone
            if field == "flex_amt_up" || field == "flex_amt_dn"
                # if input contains bus columns
                if flexible_bus_num > 0
                    # find if the flexibility of any bus is also specified in a zone column
                    bus_cols_zone_id = case.bus_zone[findall(
                        x -> x in flexible_bus_id, case.busid
                    )]
                    bus_cols_zone_idx = [
                        findfirst(isequal(i), flexible_zone_id) for i in bus_cols_zone_id
                    ]

                    # substract individual bus columns from zone total flexibility
                    for i in 1:flexible_bus_num
                        # substract from total if the zone of this bus is specified in a zone column
                        if !isnothing(bus_cols_zone_idx[i])
                            fh[:, zone_columns_idx[bus_cols_zone_idx[i]] + 1] -= fh[
                                :, bus_columns_idx[i] + 1
                            ]
                        end
                    end

                    # check if flexibility in any zone is less than the sum of buses
                    for i in 1:flexible_zone_num
                        if any(x -> x < 0, fh[:, zone_columns_idx[flexible_zone_num] + 1])
                            throw(
                                ErrorException(
                                    "Input ERROR: Total zone flexibility less than sum of bus
                       flexibility for zone " * string(flexible_zone_id[i]),
                                ),
                            )
                        end
                    end
                end

                # convert zone-level to bus-level
                converted_zone_flexibility =
                    Matrix(fh[:, [i + 1 for i in zone_columns_idx]]) * zone_to_bus_shares[zone_cols_idx, :]

                # add bus-level columns on top of converted bus-level flexibility matrix
                flex_bus_idx = zeros(Int64, 0)
                for i in sets.bus_idx
                    bus = case.busid[i]
                    # add specified flexibility to the corresponding bus
                    if bus in flexible_bus_id
                        converted_zone_flexibility[:, i] += fh[string(bus)]
                    end

                    # identify flexible bus by their total flexibility and append to list
                    if any(x -> x > 0, converted_zone_flexibility[:, i])
                        append!(flex_bus_idx, i)
                    end
                end
                eq_bus_df = DataFrames.DataFrame(
                    converted_zone_flexibility[:, flex_bus_idx]
                )
                # for cost, the zone numbers apply to all buses except those with dedicated columns
            else
                # convert zone-level cost using incidence matrix
                converted_zone_cost =
                    Matrix(fh[:, [i + 1 for i in zone_columns_idx]]) * zone_to_bus_incidence[zone_cols_idx, :]

                # replace bus-level cost for bus columns in converted bus-level cost matrix
                flex_bus_idx = zeros(Int64, 0)
                for i in sets.bus_idx
                    bus = case.busid[i]
                    # add specified cost to the corresponding bus
                    if bus in flexible_bus_id
                        converted_zone_cost[:, i] = fh[string(bus)]
                    end

                    # identify flexible bus by their total cost and append to list
                    if any(x -> x > 0, converted_zone_cost[:, i])
                        append!(flex_bus_idx, i)
                    end
                end
                eq_bus_df = DataFrames.DataFrame(converted_zone_cost[:, flex_bus_idx])
            end

            # add header and datetime index column
            DataFrames.rename!(eq_bus_df, Symbol.(case.busid[flex_bus_idx]))
            DataFrames.insertcols!(eq_bus_df, 1, :"UTC Time" => fh["UTC Time"])

            # store to new df object
            demand_flexibility_updated[field] = eq_bus_df

            # if all columns are buses, use the original dataframe
        else
            demand_flexibility_updated[field] = fh
        end
    end

    # Convert Dict to type DemandFlexibility
    demand_flexibility_updated = (;
        (Symbol(k) => v for (k, v) in demand_flexibility_updated)...
    )
    demand_flexibility_updated = DemandFlexibility(; demand_flexibility_updated...)

    return demand_flexibility_updated
end

ext: jl


text:
# Copyright (c) 2014, Dreki Þórgísl <dreki@billo.systems>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#  * Redistributions of source code must retain the above copyright notice,
#    this list of conditions and the following disclaimer.
#  * Redistributions in binary form must reproduce the above copyright notice,
#    this list of conditions and the following disclaimer in the documentation
#    and/or other materials provided with the distribution.
#  * Neither the name of the copyright holders nor the names of its
#    contributors may be used to endorse or promote products derived from this
#    software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
module Decode

using ErlPort.Exceptions

export decode, decodeterm, decodestring, decodeatom,
decodesmallint, decodeint, decodebin, decodefloat,
decodesmalltuple,
decompressterm, int2unpack, int4unpack

include("Tags.jl")
include("Util.jl")

function decode(bytes::Array{Uint8,1})
    lencheck(bytes, length(bytes) == 0)
    if bytes[1] != version
        throw(UnknownProtocolVersion(bytes[1]))
    end
    if length(bytes) >= 2 && bytes[2] == compressedtag
        # XXX maybe have this match the call to decode below? bytes[2:end]
        # instead of just bytes?
        return decodeterm(decompressterm(bytes))
    end
    return decodeterm(bytes[2:end])
end

function decode(unsupported)
    throw(UnsupportedType(unsupported))
end

function decodeterm(bytes::Array{Uint8,1})
    if length(bytes) == 0
        throw(IncompleteData(bytes))
    end
    tag = bytes[1]
    if tag == atomtag
        return decodeatom(bytes)
    elseif tag == niltag
        return decodenil(bytes)
    elseif tag == stringtag
        return decodestring(bytes)
    elseif tag == smalltupletag
        return decodesmalltuple(bytes)
    elseif tag in [listtag, smalltupletag, largetupletag]
        return bytes
    elseif tag == smallinttag
        return decodesmallint(bytes)
    elseif tag == inttag
        return decodeint(bytes)
    elseif tag == bintag
        return decodebin(bytes)
    elseif tag == newfloattag
        return decodefloat(bytes)
    elseif tag in [smallbiginttag, largebiginttag]
        # XXX move logic out to function
        if tag == smallbiginttag
            (len, sign) = (0, 0)
            tail = bytes[4:end]
        else
            (len, sign) = (0, 0)
            tail = bytes[7:end]
        end
        n = 0
        if len
            n = 0
            if sign
                n = -n
            end
        end
        return n, tail[len+1:end]
    else
        throw(UnsupportedData(bytes))
    end
end

function decodeterm(acc::Array, byte::Uint8)
    vcat(acc, decodeterm([byte]))
end

function decodeatom(bytes::Array{Uint8,1})
    len = lencheck(bytes, 3)
    unpackedlen = lencheck(len, int2unpack(bytes[2:3]) + 3, bytes)
    name = bytes[4:unpackedlen]
    if name == b"true"
        return (true, bytes[unpackedlen+1:end])
    elseif name == b"false"
        return (false, bytes[unpackedlen+1:end])
    elseif name == b"undefined"
        return (nothing, bytes[unpackedlen+1:end])
    else
        return (symbol(name), bytes[unpackedlen+1:end])
    end
end

function decodenil(bytes::Array{Uint8,1})
    #return (nothing, bytes[2:end])
    return ([], bytes[2:end])
end

function decodestring(bytes::Array{Uint8,1})
    len = lencheck(bytes, 3)
    unpackedlen = lencheck(len, int2unpack(bytes[2:3]) + 3, bytes)
    (bytes[4:unpackedlen], bytes[unpackedlen+1:end])
end

function decodesmallint(bytes::Array{Uint8,1})
    lencheck(bytes, 2)
    (bytes[2], bytes[3:end])
end

function decodeint(bytes::Array{Uint8,1})
    lencheck(bytes, 5)
    (int4unpack(bytes[2:5]), bytes[6:end])
end

function decodebin(bytes::Array{Uint8,1})
    len = lencheck(bytes, 5)
    unpackedlen = lencheck(len, int4unpack(bytes[2:5]) + 5, bytes)
    (bytes[6:unpackedlen], bytes[unpackedlen+1:end])
end

function decodefloat(bytes::Array{Uint8,1})
    lencheck(bytes, 9)
    (floatunpack(bytes[2:9]), bytes[10:end])
end

function converttoarray(len::Int64, tail::Array{Uint8,1})
    for i in len:-1:0
        (term, tail) = decodeterm(tail)
    end
end

function decodesmalltuple(bytes::Array{Uint8,1})
    lencheck(bytes, 2)
    (len, tail) = (bytes[2], bytes[3:end])
end

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using StatsBase
using Distributions
using Plots

X = randn(100)
plot(x -> ecdf(X)(x), -5, 5; label="ecdf", legend=:topleft)
plot!(x -> cdf(Normal(), x); label="Normal()")

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
@testset "Enums" begin
    test(EnumDefinition, enum_by_name, :VkAttachmentLoadOp, :(
        @cenum AttachmentLoadOp::UInt32 begin
            ATTACHMENT_LOAD_OP_LOAD = 0
            ATTACHMENT_LOAD_OP_CLEAR = 1
            ATTACHMENT_LOAD_OP_DONT_CARE = 2
        end
    ))

    test(EnumDefinition, enum_by_name, :VkResult, :(
        @cenum Result::Int32 begin
            SUCCESS = 0
            NOT_READY = 1
            TIMEOUT = 2
            EVENT_SET = 3
            EVENT_RESET = 4
            INCOMPLETE = 5
            ERROR_OUT_OF_HOST_MEMORY = -1
            ERROR_OUT_OF_DEVICE_MEMORY = -2
            ERROR_INITIALIZATION_FAILED = -3
            ERROR_DEVICE_LOST = -4
            ERROR_MEMORY_MAP_FAILED = -5
            ERROR_LAYER_NOT_PRESENT = -6
            ERROR_EXTENSION_NOT_PRESENT = -7
            ERROR_FEATURE_NOT_PRESENT = -8
            ERROR_INCOMPATIBLE_DRIVER = -9
            ERROR_TOO_MANY_OBJECTS = -10
            ERROR_FORMAT_NOT_SUPPORTED = -11
            ERROR_FRAGMENTED_POOL = -12
            ERROR_UNKNOWN = -13
            ERROR_OUT_OF_POOL_MEMORY = -1000069000
            ERROR_INVALID_EXTERNAL_HANDLE = -1000072003
            ERROR_FRAGMENTATION = -1000161000
            ERROR_INVALID_OPAQUE_CAPTURE_ADDRESS = -1000257000
            ERROR_SURFACE_LOST_KHR = -1000000000
            ERROR_NATIVE_WINDOW_IN_USE_KHR = -1000000001
            SUBOPTIMAL_KHR = 1000001003
            ERROR_OUT_OF_DATE_KHR = -1000001004
            ERROR_INCOMPATIBLE_DISPLAY_KHR = -1000003001
            ERROR_VALIDATION_FAILED_EXT = -1000011001
            ERROR_INVALID_SHADER_NV = -1000012000
            ERROR_INVALID_DRM_FORMAT_MODIFIER_PLANE_LAYOUT_EXT = -1000158000
            ERROR_NOT_PERMITTED_EXT = -1000174001
            ERROR_FULL_SCREEN_EXCLUSIVE_MODE_LOST_EXT = -1000255000
            THREAD_IDLE_KHR = 1000268000
            THREAD_DONE_KHR = 1000268001
            OPERATION_DEFERRED_KHR = 1000268002
            OPERATION_NOT_DEFERRED_KHR = 1000268003
            PIPELINE_COMPILE_REQUIRED_EXT = 1000297000
        end
    ))

    test(EnumDefinition, enum_by_name, :VkPipelineCacheHeaderVersion, :(
        @cenum PipelineCacheHeaderVersion::UInt32 begin
            PIPELINE_CACHE_HEADER_VERSION_ONE = 1
    end
    ))
end


text:
using Wasmer

using Test
using Downloads

mktempdir() do tmp_dir
    # TODO: Find a better way to provide the files
    location = joinpath(tmp_dir, "cowsay.wasm")
    Downloads.download(
        "https://registry-cdn.wapm.io/contents/_/cowsay/0.2.0/target/wasm32-wasi/release/cowsay.wasm",
        location,
    )

    wasm_code = read(location) |> Wasmer.WasmVec

    engine = WasmEngine()
    store = WasmStore(engine)
    wasm_module = WasmModule(store, wasm_code)

    msg = "Meuuuh!"

    wasi_config = Wasmer.WasiConfig("wasi_unstable")
    Wasmer.wasi_config_inherit_stderr(wasi_config)
    Wasmer.wasi_config_inherit_stdin(wasi_config)
    Wasmer.wasi_config_capture_stdout(wasi_config)

    env = Wasmer.WasiEnv(wasi_config)

    imports = Wasmer.get_imports(store, wasm_module, env)
    instance = WasmInstance(store, wasm_module, imports)

    exports(instance)._start()

    buffer = Wasmer.wasi_env_read_stdout(env, 512)
    msg = buffer .|> Char |> join |> Text

    @info "WASI worked" msg
end


text:
#SQLite

@info "SQLite.jl is available - including functionality to read / write SQLite databases"

using .SQLite

function read_table(::SQLiteFormat, filename:: AbstractString, tablename:: AbstractString; kwargs...)
    db = SQLite.DB(filename)
    return read_table(db, tablename; kwargs...)
end

"""
    read_table(db:: SQLite.DB, tablename:: AbstractString; kwargs...)

Reads the content of a whole table and returns a Tables.jl compatible object.
This method takes an instance of an SQLite database connection as input.
"""
function read_table(db:: SQLite.DB, tablename:: AbstractString; kwargs...)
     _checktablename(tablename)
    return DBInterface.execute(db, "select * from $tablename") # SQL parameters cannot be used for table names
end

read_sql(db:: SQLite.DB, sql:: AbstractString) = DBInterface.execute(db, sql)

function write_table!(::SQLiteFormat, filename:: AbstractString, tablename:: AbstractString, table; kwargs...)
    _checktable(table)
    db = SQLite.DB(filename)
    write_table!(db, tablename, table; kwargs...)
end

"""
    write_table(db:: SQLite.DB, tablename:: AbstractString, table; kwargs...)

Writes the data into a new table in the SQLite database.
As default, this operation will fail if the table already exists. Different behavior can be specified using
the `kwargs` passed to `SQLite.load!`.

Note that this behavior is different to the one for PostgreSQL, where the table must already exist.
"""
function write_table!(db:: SQLite.DB, tablename:: AbstractString, table; kwargs...)
    _checktable(table)
    table |> SQLite.load!(db, tablename; kwargs...)
    nothing
end


text:
using Images
function calc()
    width, height = 1000, 800
    img = zeros(height,width)
    a = .9
    b = -.6013
    c = 2.
    d = .8
    x = -.72
    y = -.64
    for i in 1:1000000
        x, y = x ^ 2 - y ^ 2 + a * x + b * y, 2x*y + c*x + d*y
        img[height-trunc(Int, height*(y+2)/3), trunc(Int, width*(x+1.4)/2.0)] = 1
    end
    save("tinkerbell1.png", Gray.(img))
end
@time calc()

max_stars_repo_name: BoardTM/fractals

max_issues_repo_name: BoardTM/fractals

max_forks_repo_name: BoardTM/fractals


text:
# Now test the constraint on the constraint manager

using TrajOptSOCPs
using Test

@testset "Simple Constraint Manager Base" begin

    # Make three constraints...
    c1 = AL_AffineEquality([5], [10])
    c2 = AL_AffineEquality([5], [10])
    c3 = AL_AffineEquality([5], [10])
    # ... and the associated dual
    lambdaVec = [[2], [3], [4]]

    # Make the constraint manager
    cM1 = TrajOptSOCPs.constraintManager_Base([c1, c2, c3], lambdaVec)
    x0Test = [5]
    penaltyTest = 10.0

    # Evaluate the constraints
    @test evalConstraints(cM1, x0Test, penaltyTest) == 297.0

    # Evaluate Jacobian of the Constraint Terms
    gcm = evalGradConstraints(cM1, x0Test, penaltyTest)
    @test size(gcm) == ()

end

@testset "2D Equalities and Inequalities Mixed Case" begin

    # Make three constraints of mixed types...
    c1 = AL_AffineInequality([5 0; 0 6], [10; 12])
    c2 = AL_AffineEquality([5 0; 0 6], [10; 12])
    c3 = AL_AffineInequality([5 0; 0 6], [10; 12])
    # ... and the associated dual
    lambdaVec = [[2; 2], [3; 4], [4; 2]]

    # Make the constraint manager
    cM1 = TrajOptSOCPs.constraintManager_Base([c1, c2, c3], lambdaVec)
    x0Test = [5; -3]
    penaltyTest = 10.0

    # Evaluate the constraints
    @test evalConstraints(cM1, x0Test, penaltyTest) == 567.0

    # Evaluate Gradient of the Constraint Terms
    gcm = evalGradConstraints(cM1, x0Test, penaltyTest)
    @test size(gcm, 1) == 2

    # Evaluate Hessian of the Constraint Terms
    hcm = evalHessConstraints(cM1, x0Test)
    @test size(hcm) == (2, 2)

end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
immutable NMR
  jobid::Int
  inputs::Array{String}
  fun::Function
  outputfilename::String
end

type MapperException <: Exception
  msg::String
end

type UDFException <: Exception
  msg::String
end

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ 5eebac68-891b-11eb-37ba-c1d9481c6134
begin
    using Pkg; Pkg.activate("MLJ_env", shared=true)
	using Test
	using Random
	using Distributions
	using PlutoUI
	using Printf
	using JSON

	push!(LOAD_PATH, "./src")
	using YaCounter
end

# ╔═╡ 87e5e2ec-8915-11eb-362b-6bf13a36b8e4
md"""
## NLP

ref. from book **"Data Science from Scratch"**, Chap 21

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ cb25501a-8915-11eb-3626-a1ae6d233f92
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ cb0838ae-8915-11eb-082f-3991192a101f
html"""

"""

# ╔═╡ caf03ada-8915-11eb-2935-01e110accba6
html"""

"""

# ╔═╡ 09066aa6-9049-11eb-27d6-b10b3120c900
begin
	const Int = Integer
	const T_Int = Tuple{Int, Int}
	const D_TV = Dict{Tuple{Int, Int}, Vector{Int}}
	const F = Float64
	const V_F = Vector{F}
end

# ╔═╡ cad66356-8915-11eb-0173-ada200678f63
md"""
### Gibbs Sampling

Generating samples from some distributions is easy. We can get uniform random variables with: `Uniform()`

and normal random variables with: `Normal(...)`

But some distributions are harder to sample from.

Gibbs sampling is a technique for generating samples from multidimensional distributions when we only know some of the conditional distributions.

Imagine rolling two dice. Let x be the value of the first die and y be the sum of the dice, and imagine we wanted to generate lots of (x, y) pairs. In this case it is easy to generate the samples directly:
"""

# ╔═╡ c11e47b0-9044-11eb-1f22-9de42ad879f1
Random.seed!(17)

# ╔═╡ dc3b7280-8915-11eb-258a-1789476897b9
function roll_die()::Int
	rand(1:6)
end

# ╔═╡ dc209b9a-8915-11eb-0f63-ab3b7c33793c
function direct_sample()::T_Int
	d1, d2 = roll_die(), roll_die()
	(d1, d1 + d2)
end

# ╔═╡ d7aade44-891b-11eb-00d3-abe7406c4b09
md"""
Imagine that we only knew the conditional distributions.

The distribution of y conditional on x is easy — if we know the value of x, y is
equally likely to be x + 1, x + 2, x + 3, x + 4, x + 5, or x + 6:
"""

# ╔═╡ 9ae99efe-891c-11eb-3d1f-89182e3f9ff8
function random_y_given_x(x::Integer)::Int
	"""
	Equally likely to be x + 1, x + 2, ... , x + 6
	"""
	x + roll_die()
end

# ╔═╡ ea1f4a2c-891d-11eb-3027-bd447dd9d1e1
md"""
However the other direction is a bit more complicated. For example, if we know that y is
2, then necessarily x is 1 (since the only way two dice can sum to 2 is if
both of them are 1). If we know y is 3, then x is equally likely to be 1 or 2.
Similarly, if y is 11, then x has to be either 5 or 6:
"""

# ╔═╡ 13e3ccb6-8aae-11eb-3a9e-85b5afc74097
function random_x_given_y(y::Integer)::Int
	if y ≤ 7
		## if total is less (or equal) than 7, the first die is equally likely to be
		## 1, 2, ..., (total - 1)
		rand(1:y)
	else
		## if the total is 7+, the first die is equally likely to be
		## (total - 6), (total - 5), ..., 6
		rand(y - 6:7)
	end
end

# ╔═╡ 13c729da-8aae-11eb-29aa-b171020ad8dd
md"""
The way Gibbs sampling works is that we start with any (valid) values for x and y and then repeatedly alternate:
  - replacing x with a random value picked conditional on y and
  - replacing y with a random value picked conditional on x.

After a number of iterations, the resulting values of x and y will represent a sample from the unconditional joint distribution:
"""

# ╔═╡ 2e9609b2-9046-11eb-0ac6-7ff568ff8668
function gibbs_sample(;num_iters=100)::T_Int
	x, y = 1, 2 # doesn't really matter
	for _ ∈ 1:num_iters
		x = random_x_given_y(y)
		y = random_y_given_x(x)
	end
	(x, y)
end

# ╔═╡ 13b2940c-8aae-11eb-2dd2-8dec30dd4881
function cmp_distributions(;num_samples=1000)::D_TV
	counts = D_TV()
	for _ ∈ 1:num_samples
		gs = gibbs_sample()
		v_s = get(counts, gs, [0, 0])
		v_s[1] += 1
		counts[gs] = v_s
		ds = direct_sample()
		v_s = get(counts, ds, [0, 0])
		v_s[2] += 1
		counts[ds] = v_s
	end
	counts
end

# ╔═╡ 8aa8bb6a-9047-11eb-1c48-7dfd972d7580
cmp_distributions()

# ╔═╡ dc0c5a5e-8915-11eb-2c62-5328b42add69
html"""

"""

# ╔═╡ dbefe77a-8915-11eb-39f1-5d0c84729739
md"""
### Topic Modeling

An approach to understanding our users’ interests might try to identify the topics
that underlie those interests. A technique called latent Dirichlet allocation (LDA) is commonly used
to identify common topics in a set of documents.
We’ll apply it to documents that consist of each user’s interests.

LDA has some similarities to the Naive Bayes classifier we built previously in that it assumes a probabilistic
model for documents. For our purposes the model assumes that:

  - There is some fixed number K of topics.

  - There is a random variable that assigns each topic an associated probability distribution over words.
    Think of this distribution as the probability of seeing word w given topic k.

  - There is another random variable that assigns each document a probability distribution over topics.
    Think of this distribution as the mixture of topics in document d.

  - Each word in a document was generated by first randomly picking a topic (from the document’s distribution
    of topics) and then randomly picking a word (from the topic’s distribution of words).

In particular, we have a collection of documents, each of which is a list of words. And we have a corresponding
collection of document_topics that assigns a topic (here a number between 0 and K – 1) to each word in each
document.

So, the fifth word in the fourth document is: documents[4][5]
and the topic from which that word was chosen is: document_topics[4][5]

This very explicitly defines each document’s distribution over topics, and it implicitly defines each topic’s distribution over words.

We can estimate the likelihood that topic 1 produces a certain word by comparing how many times topic 1 produces that word with how many
times topic 1 produces any word. (Similarly, when we built a spam filter, we compared how many times each word appeared in spams with the
total number of words appearing in spams.)

Although these topics are just numbers, we can give them descriptive names by looking at the words on which they put the heaviest weight.
We just have to somehow generate the document_topics. This is where Gibbs sampling comes into play.

  We start by assigning every word in every document a topic completely at random.

  Now we go through each document one word at a time. For that word and document, we construct weights for each topic that depend on the
  (current) distribution of topics in that document and the (current) distribution of words for that topic.

  We then use those weights to sample a new topic for that word.

  If we iterate this process many times, we will end up with a joint sample from the topic-word distribution and the document-topic distribution.

Let us start with a function to randomly choose an index based on an arbitrary set of weights:
"""

# ╔═╡ dbd86924-8915-11eb-351f-8362f09ba984
function sample_from(weights::V_F)::Int
	"""
	returns ix (index) with probability weights[ix] / sum(weights)
	"""
	total = sum(weights)
	rnd = total * rand()       ## uniform between 0 and total
	for (ix, w) ∈ enumerate(weights)
		rnd -= w
		rnd ≤ 0. && (return ix) ## return the smallest i such that
		#                       ## weights[0] + ... + weights[i] ≥ rnd
	end
end

# ╔═╡ d32b86e0-8928-11eb-193c-e3d7c85d23dd
md"""
if we give it weights [1, 1, 3] , then one-fifth of the time it will return 0, one-fifth of the time it will return 1, and three-fifths of the time it will return 2.
"""

# ╔═╡ 4a92cc8e-904a-11eb-11de-e77d09365187
begin
	draws = Counter([sample_from([.1, .1, .8]) for _ ∈ 1:1000])

	@test 10 < draws[1] < 190    ## should be ~10%, this is a really loose test
	@test 10 < draws[2] < 190    ## should be ~10%, this is a really loose test
	@test 650 < draws[3] < 950   ## should be ~80%, this is a really loose test
	@test sum(values(draws)) == 1000
end

# ╔═╡ f1bc0f20-904a-11eb-394e-23d1ee736caf
begin
	const Documents = [
    	["Hadoop", "Big Data", "HBase", "Java", "Spark", "Storm", "Cassandra"],
    	["NoSQL", "MongoDB", "Cassandra", "HBase", "Postgres"],
    	["Python", "scikit-learn", "scipy", "Julia",  "numpy", "statsmodels", "pandas"],
    	["R", "Python", "Julia", "statistics", "regression", "probability"],
    	["machine learning", "regression", "decision trees", "libsvm"],
    	["Python", "R", "Java", "C++", "Julia", "programming languages"],
    	["statistics", "probability", "mathematics", "theory"],
    	["machine learning", "scikit-learn", "Mahout", "neural networks"],
    	["neural networks", "deep learning", "Flux", "Big Data", "artificial intelligence"],
    	["Hadoop", "Java", "MapReduce", "Big Data"],
    	["statistics", "R", "statsmodels"],
    	["C++", "deep learning", "artificial intelligence", "probability"],
    	["pandas", "DataFrames", "Julia", "R", "Python"],
    	["databases", "HBase", "Postgres", "MySQL", "MongoDB"],
    	["libsvm", "regression", "support vector machines"],
		["Julia", "MLJ", "Flux", "DataFrames"]
	]

	const K = 4;
end

# ╔═╡ f1a241d0-904a-11eb-06cc-fda025d4f3b5
begin
	## a list of Counters, one for each document
	doc_topic_counts = [Counter() for _ ∈ Documents]

	## a list of Counters, one for each  topic
	topic_word_counts = [Counter() for _ ∈ 1:K]

	topic_counts = zeros(Int, K)     # a vector of numbers, one for each topic
	doc_lengths = length.(Documents) # a vector of numbers, one for each document

	unique_words = [w for doc ∈ Documents for w ∈ doc] |> unique
	const W, D = length(unique_words), length(Documents);
end

# ╔═╡ d1f14cce-904e-11eb-27c9-4353466db275
md"""
Now we are ready to define our conditional probability functions. Each Will have  smoothing term that ensures every topic has a nonzero chance of being chosen in any document and that every word has a nonzero chance of being chosen for any topic:
"""

# ╔═╡ d1b8b198-904e-11eb-0a40-27d960d3b303
function p_topic_given_doc(topic::Int, d::Int; ϵ=.1)::F
	"""
	A closure calculating the fraction of words in document 'd' that are assigned
	to 'topic' (plus some smoothing)
	"""
	(doc_topic_counts[d][topic] + ϵ) / (doc_lengths[d] + K * ϵ)
end

# ╔═╡ 60b8d5ca-905d-11eb-1238-43d2ed28bdfc
function p_word_given_topic(word::String, topic::Int; ϵ=.1)::F
	"""
	The fraction of words assigned to 'topic' that equal 'word' (plus some smoothing)
	"""
	(topic_word_counts[topic][word] + ϵ) / (topic_counts[topic] + W * ϵ)
end

# ╔═╡ 4ae5db5e-904f-11eb-2073-b17993b09f0f
function topic_weight(d::Int, word::String, k::Int)::F
	"""
	Given a document and a word in that document, return the weight for the kth topic
	"""
	p_word_given_topic(word, k) * p_topic_given_doc(k, d)
end

# ╔═╡ ea0dcb42-905e-11eb-25df-05704be5e114
function choose_new_topic(d::Int, word::String)::Int
	sample_from([topic_weight(d, word, k) for k ∈ 1:K])
end

# ╔═╡ 584b4fbc-905f-11eb-0631-ad00a20dc881
md"""
Given a word and its document, the likelihood of any topic choice depends on both how likely that topic is for the document and how likely that word is for the topic.

Now let us start assigning every word to a random topic and populating our counters appropriately:
"""

# ╔═╡ 4aca19be-904f-11eb-26af-bd9f21c61a0c
begin
	Random.seed!(42)

	doc_topics = [[rand(1:K) for w ∈ doc] for doc ∈ Documents]

	for d ∈ 1:D
		for (w, topic) ∈ zip(Documents[d], doc_topics[d])
			doc_topic_counts[d][topic] += 1
			topic_word_counts[topic][w] += 1
			topic_counts[topic] += 1
		end
	end
end

# ╔═╡ 4ed6af9a-905e-11eb-1b8a-152f021dabf4
md"""
The goal is to get a joint sample of the topics-word distribution and the documents-topic distribution. We do this using a form of Gibbs sampling that uses the conditional probabilities defined previously:
"""

# ╔═╡ 4a928fa8-904f-11eb-2f47-4d319a1a6e58
for iter ∈ 1:1000
	for d ∈ 1:D, (ix, (word, topic)) ∈ enumerate(zip(Documents[d], doc_topics[d]))
		## Remove this word-topic from the counts so that it does not
		## influence the weigths
		doc_topic_counts[d][topic] -= 1
		topic_word_counts[topic][word] -= 1
		topic_counts[topic] -= 1
		doc_lengths[d] -= 1

		## Choose new topic based on the weigths
		new_topic = choose_new_topic(d, word)
		doc_topics[d][ix] = new_topic

		## Update new topic
		doc_topic_counts[d][new_topic] += 1
		topic_word_counts[new_topic][word] += 1
		topic_counts[new_topic] += 1
		doc_lengths[d] += 1
	end
end

# ╔═╡ df809854-9071-11eb-0bed-1f9e8c66b023
md"""
The topics are just numbers (from 1 to 4), if we want names we have to find by ourselves.
"""

# ╔═╡ f3d142e2-906f-11eb-1a04-37c8d25ac050
with_terminal() do
	# 7 items per topic
	for (k, word_cnts) ∈ enumerate(topic_word_counts)
		for (ix, (w, cnt)) ∈ enumerate(most_common(word_cnts, 7))
			cnt > 0 && @printf("%3d/%2d - %-25s %4d\n", k, ix, w, cnt)
		end
		println()
	end
end

# ╔═╡ 4f708e46-9070-11eb-2e90-45e2e028bbdb
md"""
We may define the following topics:
  1. Big Data & Databases
  2. AI and Languages
  3. ML
  4. Statistics and Dynamic Languages
"""

# ╔═╡ 4b449308-9070-11eb-2aed-65775ffabadb
topic_names = [
	"Big Data & Databases",
	"AI and Languages",
	"ML",
	"Statistics and Dynamic Languages"
]

# ╔═╡ 89f83412-9073-11eb-11bf-e7a77ea377ab
with_terminal() do
	for (doc, topic_cnts) ∈ zip(Documents, doc_topic_counts)
		println(doc)
		for (topic, cnt) ∈ most_common(topic_cnts)
			cnt > 0. && @printf(" - %-20s => %2d\n", topic_names[topic], cnt)
		end
		println()
	end
end

# ╔═╡ Cell order:
# ╟─87e5e2ec-8915-11eb-362b-6bf13a36b8e4
# ╠═5eebac68-891b-11eb-37ba-c1d9481c6134
# ╟─cb25501a-8915-11eb-3626-a1ae6d233f92
# ╟─cb0838ae-8915-11eb-082f-3991192a101f
# ╟─caf03ada-8915-11eb-2935-01e110accba6
# ╠═09066aa6-9049-11eb-27d6-b10b3120c900
# ╟─cad66356-8915-11eb-0173-ada200678f63
# ╠═c11e47b0-9044-11eb-1f22-9de42ad879f1
# ╠═dc3b7280-8915-11eb-258a-1789476897b9
# ╠═dc209b9a-8915-11eb-0f63-ab3b7c33793c
# ╟─d7aade44-891b-11eb-00d3-abe7406c4b09
# ╠═9ae99efe-891c-11eb-3d1f-89182e3f9ff8
# ╟─ea1f4a2c-891d-11eb-3027-bd447dd9d1e1
# ╠═13e3ccb6-8aae-11eb-3a9e-85b5afc74097
# ╟─13c729da-8aae-11eb-29aa-b171020ad8dd
# ╠═2e9609b2-9046-11eb-0ac6-7ff568ff8668
# ╠═13b2940c-8aae-11eb-2dd2-8dec30dd4881
# ╠═8aa8bb6a-9047-11eb-1c48-7dfd972d7580
# ╟─dc0c5a5e-8915-11eb-2c62-5328b42add69
# ╟─dbefe77a-8915-11eb-39f1-5d0c84729739
# ╠═dbd86924-8915-11eb-351f-8362f09ba984
# ╟─d32b86e0-8928-11eb-193c-e3d7c85d23dd
# ╠═4a92cc8e-904a-11eb-11de-e77d09365187
# ╠═f1bc0f20-904a-11eb-394e-23d1ee736caf
# ╠═f1a241d0-904a-11eb-06cc-fda025d4f3b5
# ╟─d1f14cce-904e-11eb-27c9-4353466db275
# ╠═d1b8b198-904e-11eb-0a40-27d960d3b303
# ╠═60b8d5ca-905d-11eb-1238-43d2ed28bdfc
# ╠═4ae5db5e-904f-11eb-2073-b17993b09f0f
# ╠═ea0dcb42-905e-11eb-25df-05704be5e114
# ╟─584b4fbc-905f-11eb-0631-ad00a20dc881
# ╠═4aca19be-904f-11eb-26af-bd9f21c61a0c
# ╟─4ed6af9a-905e-11eb-1b8a-152f021dabf4
# ╠═4a928fa8-904f-11eb-2f47-4d319a1a6e58
# ╟─df809854-9071-11eb-0bed-1f9e8c66b023
# ╠═f3d142e2-906f-11eb-1a04-37c8d25ac050
# ╟─4f708e46-9070-11eb-2e90-45e2e028bbdb
# ╠═4b449308-9070-11eb-2aed-65775ffabadb
# ╠═89f83412-9073-11eb-11bf-e7a77ea377ab

max_stars_repo_path: Data Science From Scratch/21-NLP.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Data Science From Scratch/21-NLP.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Data Science From Scratch/21-NLP.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
function SelectVariable(VariableDef,X_Tilde, delta1 ,delta2, MU ;
                  n1 = [], n2 = [], TAV = 0,ZZ = 0,xlb = [],xub = [],Type=1)

  if Type == 1 || isempty(Type) # Randomly Chosen
    NonIntegers = find(.!isinteger.(X_Tilde[VariableDef[1]+1:end]))
    return VariableDef[1]+NonIntegers[rand(1:length(NonIntegers))]
  elseif Type == 2 # Randomly among those which have the most value violation.
    NonIntegers = find(.!isinteger.(X_Tilde[VariableDef[1]+1:end]))
    MaxInd = findmax( min.(X_Tilde[VariableDef[1]+NonIntegers]-floor.(X_Tilde[VariableDef[1]+NonIntegers]),
                      ceil.(X_Tilde[VariableDef[1]+NonIntegers])-X_Tilde[VariableDef[1]+NonIntegers]) )[2];
    return VariableDef[1] + MaxInd ;
  elseif Type == 3
    NonIntegers = find(.!isinteger.(X_Tilde[VariableDef[1]+1:end]))
    S = MU*min.(delta1,delta2) + (1-MU)*max.(delta1,delta2)

    ind = findmax(S[NonIntegers])[2]
    # println(ind)
    return NonIntegers[ind]+VariableDef[1]

  elseif Type == 4
    NonIntegers = find(.!isinteger.(X_Tilde[VariableDef[1]+1:end]))
    S = zeros(length(NonIntegers))
    for i = 1:length(NonIntegers)
      NewXUB = copy(xub);
      NewXUB[VariableDef[1]+NonIntegers[i]] = floor(X_Tilde[VariableDef[1]+NonIntegers[i] ])
      Z2 = RSolver(A,b,C,d,xlb,NewXUB)[1];

      NewXLB = copy(xlb);
      NewXLB[VariableDef[1]+NonIntegers[i]] = ceil(X_Tilde[VariableDef[1]+NonIntegers[i]])
      Z1 = RSolver(A,b,C,d,NewXLB,xub)[1];

      if n1[NonIntegers[i]] < TAV && !isempty(Z1)
        Delta1 = Z1 - ZZ;
      else
        Delta1 = delta1[i]
      end
      if n2[NonIntegers[i]] < TAV && !isempty(Z2)
        Delta2 = Z2 - ZZ;
      else
        Delta2 = delta2[i]
      end
      S[i] = MU*min(Delta1, Delta2) + (1-MU)*max(Delta1,Delta2)
      ind = findmax(S)[2]

      return NonIntegers[ind]+VariableDef[1]

    end
  end

end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.DLR_Visualization (Modia3D/renderer/DLR_Visualization/_module.jl)
#

### Functions available for community and professional edition of SimVis

function SimVis_init(simVis::SimVis_Renderer)
   ccall(simVis.init, Nothing,(Cstring,Cstring,Cint,Cint),
                               simVis.directory, simVis.host, simVis.port, Int(simVis.sync))
end

function SimVis_shutdown(simVis::SimVis_Renderer)
   ccall(simVis.shutdown, Nothing,())
end

function SimVis_getObjectID(simVis::SimVis_Renderer, emptyObjectID::Int)
   ccall(simVis.getObjectID, Ptr{Nothing},(Cint,), emptyObjectID)
end

function SimVis_freeObjectID(simVis::SimVis_Renderer, obj::Ptr{Nothing})
   ccall(simVis.freeObjectID, Nothing,(Ptr{Nothing},), obj)
end

function SimVis_setTime(simVis::SimVis_Renderer, time::Float64)
   ccall(simVis.setTime, Nothing,(Cdouble,), time)
end

function SimVis_setBaseObject(simVis::SimVis_Renderer,
                              ID::Ptr{Nothing},
                              state::Cint,
                              baseObjType::Cint,
                              pos::MVector{3,Float64},
                              T::MMatrix{3,3,Float64,9},
                              scale::MVector{3,Float64},
                              color::MVector{3,Cint},
                              wireframe::Cint,
                              reflectslight::Cint,
                              specularCoefficient::Float64,
                              extra::MVector{3,Float64},
                              alpha::Float64,
                              canCollide::Cint,
                              shadowMask::Cint)
   ccall(simVis.setBaseObject, Nothing,
           (Ptr{Nothing},Cint,Cint,Ref{SVector{3,Float64}},Ref{SMatrix{3,3,Float64,9}},Ref{MVector{3,Float64}},Ref{MVector{3,Cint}},Cint,Cint,Cdouble,Ref{MVector{3,Float64}},Cdouble,Cint,Cint),
           ID,state,baseObjType,pos,T,scale,color,wireframe,reflectslight,specularCoefficient,extra,alpha,canCollide,shadowMask)
end

function SimVis_setBaseObject(simVis::SimVis_Renderer,
                              ID::Ptr{Nothing},
                              state::Cint,
                              baseObjType::Cint,
                              pos::SVector{3,Float64},
                              T::SMatrix{3,3,Float64,9},
                              scale::MVector{3,Float64},
                              color::MVector{3,Cint},
                              wireframe::Cint,
                              reflectslight::Cint,
                              specularCoefficient::Float64,
                              extra::MVector{3,Float64},
                              alpha::Float64,
                              canCollide::Cint,
                              shadowMask::Cint)
   ccall(simVis.setBaseObject, Nothing,
           (Ptr{Nothing},Cint,Cint,Ref{SVector{3,Float64}},Ref{SMatrix{3,3,Float64,9}},Ref{MVector{3,Float64}},Ref{MVector{3,Cint}},Cint,Cint,Cdouble,Ref{MVector{3,Float64}},Cdouble,Cint,Cint),
           ID,state,baseObjType,pos,T,scale,color,wireframe,reflectslight,specularCoefficient,extra,alpha,canCollide,shadowMask)
end

function SimVis_setFileObject(simVis::SimVis_Renderer,
                              ID::Ptr{Nothing},
                              state::Cint,
                              pos::MVector{3,Float64},
                              T::MMatrix{3,3,Float64,9},
                              scale::MVector{3,Float64},
                              reflectslight::Cint,
                              specularCoefficient::Float64,
                              alpha::Float64,
                              wireframe::Cint,
                              canCollide::Cint,
                              filename::String,
                              smooth::Cint,
                              overwriteColor::Bool,
                              color::MVector{3,Cint},
                              shadowMask::Cint,
                              shaderName::String)
   ccall(simVis.setFileObject, Nothing,
            (Ptr{Nothing},Cint,Ref{MVector{3,Float64}},Ref{MVector{3,Float64}},Ref{MVector{3,Float64}},Cint,Cdouble,Cdouble,Cint,Cint,Cstring,Cint,Cint,Ref{MVector{3,Cint}},Cint,Cstring),
             ID,state,pos,T,scale,reflectslight,specularCoefficient,alpha,wireframe,canCollide,filename,smooth,overwriteColor,color,shadowMask,shaderName)
end

function SimVis_setFileObject(simVis::SimVis_Renderer,
                              ID::Ptr{Nothing},
                              state::Cint,
                              pos::SVector{3,Float64},
                              T::SMatrix{3,3,Float64,9},
                              scale::MVector{3,Float64},
                              reflectslight::Cint,
                              specularCoefficient::Float64,
                              alpha::Float64,
                              wireframe::Cint,
                              canCollide::Cint,
                              filename::String,
                              smooth::Cint,
                              overwriteColor::Bool,
                              color::MVector{3,Cint},
                              shadowMask::Cint,
                              shaderName::String)
   ccall(simVis.setFileObject, Nothing,
            (Ptr{Nothing},Cint,Ref{SVector{3,Cdouble}},Ref{SMatrix{3,3,Cdouble,9}},Ref{MVector{3,Float64}},Cint,Cdouble,Cdouble,Cint,Cint,Cstring,Cint,Cint,Ref{MVector{3,Cint}},Cint,Cstring),
             ID,state,pos,T, scale,reflectslight,specularCoefficient,alpha,wireframe,canCollide,filename,smooth,overwriteColor,color,shadowMask,shaderName)
end

### Functions only available for professional edition of SimVis
function SimVis_setTextObject(simVis::SimVis_Renderer,
                              ID::Ptr{Nothing},
                              screenAlignment::Cint,
                              text::String,
                              textvalue::Float64,
                              valueactive::Cint,
                              pos::MVector{3,Float64},
                              T::MMatrix{3,3,Float64,9},
                              charsize::Float64,
                              fontname::String,
                              color::MVector{3,Cint},
                              alpha::Float64,
                              offset::MVector{3,Float64},
                              alignment::Cint,
                              digits::Cint)
   ccall(simVis.setTextObject, Nothing,
           (Ptr{Nothing},Cint,Cstring,Cdouble,Cint,Ref{MVector{3,Float64}},Ref{MMatrix{3,3,Float64,9}},Cdouble,
            Cstring,Ref{MVector{3,Cint}},Cdouble,Ref{MVector{3,Float64}},Cint,Cint),
           ID, screenAlignment, text, textvalue, valueactive, pos, T, charsize,
           fontname, color, alpha, offset, alignment, digits)
end

function SimVis_setTextObject(simVis::SimVis_Renderer,
                              ID::Ptr{Nothing},
                              screenAlignment::Cint,
                              text::String,
                              textvalue::Float64,
                              valueactive::Cint,
                              pos::SVector{3,Float64},
                              T::SMatrix{3,3,Float64,9},
                              charsize::Float64,
                              fontname::String,
                              color::MVector{3,Cint},
                              alpha::Float64,
                              offset::MVector{3,Float64},
                              alignment::Cint,
                              digits::Cint)
   ccall(simVis.setTextObject, Nothing,
           (Ptr{Nothing},Cint,Cstring,Cdouble,Cint,Ref{SVector{3,Cdouble}},Ref{SMatrix{3,3,Cdouble,9}},Cdouble,
            Cstring,Ref{MVector{3,Cint}},Cdouble,Ref{MVector{3,Float64}},Cint,Cint),
           ID, screenAlignment, text, textvalue, valueactive, pos, T, charsize,
           fontname, color, alpha, offset, alignment, digits)
end


text:
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
default_background!(vis)

# box
dir = "/home/taylor/Research/ContactControl.jl/src/dynamics/quadruped/payload_mesh/Box.obj"
obj = MeshFileObject(dir)
setobject!(vis[:box], obj)
tr = Translation(0.495, 0.0, 0.3)
rt = LinearMap(0.2 * RotZ(-0.5 * π) * RotX(1.0 * π))
tform = compose(tr, rt)
settransform!(vis[:box], tform)

include_model("simple_manipulator")

# arm 1
l1 = Cylinder(Point3f0(0, 0, 0),
	Point3f0(0, 0, 1.0),
	convert(Float32, 0.025))
setobject!(vis["l1_1"], l1, MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
l2 = Cylinder(Point3f0(0, 0, 0), Point3f0(0, 0, 1.0),
	convert(Float32, 0.025))
setobject!(vis["l2_1"], l2, MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

setobject!(vis["elbow_1"], Sphere(Point3f0(0),
	convert(Float32, 0.05)),
	MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
setobject!(vis["ee_1"], Sphere(Point3f0(0),
	convert(Float32, 0.05)),
	MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))

p1 = [-0.99; -0.01; 0.0]
q1 = [0.41 * π; -0.635 * π]
p_mid = [kinematics_mid(model, q1)[1], 0.0, kinematics_mid(model, q1)[2]] + p1
p_ee = [kinematics_ee(model, q1)[1], 0.0, kinematics_ee(model, q1)[2]] + p1

settransform!(vis["l1_1"], cable_transform(p1, p_mid))
settransform!(vis["l2_1"], cable_transform(p_mid, p_ee))

settransform!(vis["elbow_1"], Translation(p_mid))
settransform!(vis["ee_1"], Translation(p_ee))

# arm 2
setobject!(vis["l1_2"], l1, MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
setobject!(vis["l2_2"], l2, MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

setobject!(vis["elbow_2"], Sphere(Point3f0(0),
	convert(Float32, 0.05)),
	MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
setobject!(vis["ee_2"], Sphere(Point3f0(0),
	convert(Float32, 0.05)),
	MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))

p2 = [0.99; -0.01; 0.0]
q2 = [0.4 *  π; 0.825 * π]
p_mid = [kinematics_mid(model, q2)[1], 0.0, kinematics_mid(model, q2)[2]] + p2
p_ee = [kinematics_ee(model, q2)[1], 0.0, kinematics_ee(model, q2)[2]] + p2

settransform!(vis["l1_2"], cable_transform(p2, p_mid))
settransform!(vis["l2_2"], cable_transform(p_mid, p_ee))

settransform!(vis["elbow_2"], Translation(p_mid))
settransform!(vis["ee_2"], Translation(p_ee))

settransform!(vis["/Cameras/default"],
	   compose(Translation(0.0, -95.0, -1.0), LinearMap(RotY(0.0 * π) * RotZ(-π / 2.0))))
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 60)

setobject!(vis["box1"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0),
		Vec(0.25, 0.5, 0.25)), MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, 1.0)))
	settransform!(vis["box1"], Translation(0.51, -0.25, 0))

setobject!(vis["box2"], GeometryBasics.HyperRectangle(Vec(0.0, 0.0, 0.0),
		Vec(0.25, 0.5, 0.25)), MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, 1.0)))
	settransform!(vis["box2"], Translation(-0.16, -0.25, 0))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
module Burgers

using NeuralOperators
using Flux
using CUDA

include("data.jl")

__init__() = register_burgers()

function train()
    if has_cuda()
        @info "CUDA is on"
        device = gpu
        CUDA.allowscalar(false)
    else
        device = cpu
    end

    modes = (16, )
    ch = 64 => 64
    σ = gelu
    m = Chain(
        Dense(2, 64),
        FourierOperator(ch, modes, σ),
        FourierOperator(ch, modes, σ),
        FourierOperator(ch, modes, σ),
        FourierOperator(ch, modes),
        Dense(64, 128, σ),
        Dense(128, 1),
        flatten
    ) |> device

    loss(𝐱, 𝐲) = sum(abs2, 𝐲 .- m(𝐱)) / size(𝐱)[end]

    loader_train, loader_test = get_dataloader()

    function validate()
        validation_losses = [loss(device(𝐱), device(𝐲)) for (𝐱, 𝐲) in loader_test]
        @info "loss: $(sum(validation_losses)/length(loader_test))"
    end

    data = [(𝐱, 𝐲) for (𝐱, 𝐲) in loader_train] |> device
    opt = Flux.Optimiser(WeightDecay(1f-4), Flux.ADAM(1f-3))
    call_back = Flux.throttle(validate, 5, leading=false, trailing=true)
    Flux.@epochs 500 @time(Flux.train!(loss, params(m), data, opt, cb=call_back))
end

end


text:
#i64o-utilities.jl

#implement the deep copy constructor and the overwriting copy! function for ArrayNums.
Base.copy{FSS}(a::ArrayNum{FSS}) = ArrayNum{FSS}(copy(a.a))
function Base.copy!{FSS}(dest::ArrayNum{FSS}, src::ArrayNum{FSS})
  for idx = 1:__cell_length(FSS)
    @inbounds dest.a[idx] = src.a[idx]
  end
end

#bits function for hlayer output.
Base.bits{FSS}(a::ArrayNum{FSS}) = mapreduce(bits, (s1, s2) -> string(s1, s2), "", a.a)

#forwarding getindex and setindex!
Base.getindex{FSS}(a::ArrayNum{FSS}, key...) = getindex(a.a, key...)
Base.setindex!{FSS}(a::ArrayNum{FSS}, X, keys...) = setindex!(a.a, X, keys...)

doc"""
`Unums.set_bit!` sets a bit in the ArrayNum referred to by the value b, this bit
is one-indexed with the bit 1 being the most significant.  A value of zero has
undefined effects.  Useful for setting bits after shifting a non-subnormal value.
"""
function set_bit!{FSS}(a::ArrayNum{FSS}, bit::UInt16)
  a_index = ((bit - o16) >> 6) + o16
  b_index = ((bit - o16) % 64)
  @inbounds a.a[a_index] = a.a[a_index] | (0x8000_0000_0000_0000 >> b_index)
  a
end

doc"""
`Unums.set_bit` sets a bit in a UInt64, this bit is one-indexed with the top bit
being the the most significant.  A value of zero has undefined effects.  Useful
for setting bits after shifting a non-subnormal value.
"""
function set_bit(a::UInt64, bit::UInt16)
  a | (0x8000_0000_0000_0000 >> (bit - o16))
end
doc"""
`Unums.frac_set_bit!(x, bit)` sets (one-indexed) bit, which is useful for setting
bits after shifting a non-subnormal value.
"""
@fracproc set_bit bit

doc"""
`Unums.get_bit(x, bit)` returns true if the (zero-indexed) bit of x is one, false
if not.
"""
function get_bit(a::UInt64, bit::UInt16)
  (a & (t64 >> bit)) != 0
end
function get_bit{FSS}(a::ArrayNum{FSS}, bit::UInt16)
  a_index = ((bit - o16) >> 6) + o16
  b_index = ((bit - o16) % 64)
  @inbounds return (a.a[a_index] & (t64 >> b_index)) != 0
end

doc"""
`Unums.copy_top(x, val)` performs the logical or of the value with the fraction
of x, or the first element in x fraction array.
"""
copy_top(a::UInt64, pattern::UInt64) = a | pattern
copy_top!{FSS}(a::ArrayNum{FSS}, pattern::UInt64) = (a[1] |= pattern)

doc"""
`Unums.frac_copy_top!(x::Unum, pattern::UInt64)` performs the logical or of the pattern
with the fraction of x, or the first element in x fraction array.
"""
@fracproc copy_top pattern

doc"""
  `Unums.inward_ubit_crosses_zero(a::UInt64, fsize)`
  `Unums.inward_ubit_crosses_zero(a::ArrayNum, fsize)`

  tells if an inward ubit from this number would cross zero (into the next exponential range)
"""
function inward_ubit_crosses_zero(a::UInt64, fsize::UInt16)
  #the mask is generated using this formula:
  mask = z64 - (o64 << (0x40 - fsize - o16))
  #check to make sure that everything is zeros
  return ((a & mask) == z64)
end
function inward_ubit_crosses_zero{FSS}(a::ArrayNum{FSS}, fsize::UInt16)
  #compute the last cell we need to scan.
  middle_spot = div(fsize, 0x0040) + 1
  middle_mask = fsize % 0x0040
  for idx = 1:(middle_spot - 1)
    @inbounds (a.a[idx] != z64) && return false
  end
  @inbounds return inward_ubit_crosses_zero(a.a[middle_spot], middle_mask)
end

#__minimum_data_width
#calculates the minimum data width to represent the passed superint.
function __minimum_data_width{FSS}(n::ArrayNum{FSS})
  res = max(z16, max_fsize(FSS) - ctz(n))
  res == 0xFFFF ? z16 : res
end
  #explanation of formula:
  #length(a) << 6:            total bits in the array representation
  #-trailing_zeros(f):        how many zeros are at the end, we can trim those
  #-1:                        the bit representation (1000...0000) = "1" has
  #                           width 0 as per our definition.
  #max(0, ...):               bit representation of (0000...0000) = "0" also
  #                           has width 0, not width "-1".

#this is a better formula for a single-width unsigned integer representation.
__minimum_data_width(n::UInt64) = (res = max(z16, 0x003F - ctz(n)); res == 0xFFFF ? z16 : res)

#simply assign this to a hash of the array itself.
Base.hash{FSS}(n::ArrayNum{FSS}, h::UInt) = hash(n.a, h)

function contract_outer_fsize(n::UInt64, s::UInt16)
  (0x0040 - ctz(~n & mask_top(s - o16))) * (s != 0)
end

function contract_inner_fsize(n::UInt64, s::UInt16)
  (0x0040 - ctz(n & mask_top(s - o16))) * (s != 0)
end

function contract_outer_fsize{FSS}(n::ArrayNum{FSS}, s::UInt16)
  #this function is dedicated to the Hot Baked Goods.
  middle_cell = div(s, 0x0040) + o16
  middle_size = s % 0x0040

  local contracted_size::UInt16
  #do our middle cell thing.
  @inbounds contracted_size = contract_outer_fsize(n.a[middle_cell], middle_size)

  for idx = middle_cell-1 : -1 : 1
    if contracted_size != 0
      return (idx * 0x0040 + contracted_size)
    else
      @inbounds contracted_size = contract_outer_fsize(n.a[idx], 0x003F)
    end
  end
  return contracted_size
end

function contract_inner_fsize{FSS}(n::ArrayNum{FSS}, s::UInt16)
  #this function is dedicated to the Hot Baked Goods.
  middle_cell = div(s, 0x0040) + o16
  middle_size = s % 0x0040

  local contracted_size::UInt16
  #do our middle cell thing.
  @inbounds contracted_size = contract_inner_fsize(n.a[middle_cell], middle_size)

  for idx = middle_cell-1 : -1 : 1
    if contracted_size != 0
      return (idx * 0x0040 + contracted_size)
    else
      @inbounds contracted_size = contract_inner_fsize(n.a[idx], 0x003F)
    end
  end
  return contracted_size
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl