text:
using DelimitedFiles, Glob, Plots

## Helper functions
#= Input demand matrix computed for all combinations of i/j and the shortest path computed,
    then return the mean load for the load blocks computed =#

##  Form a sequence containing the average demand in computed load blocks
#   Input: Mean Demand Matrix D; path= sequence of indices corresponding to optimal LB starts

function meanLoad(D, path)
    avgedLoad = []
    for i = 1:length(path)-1
        s = path[i]
        t = path[i+1]
        push!(avgedLoad, D[s, t])
    end
    return avgedLoad
end

##  Shortest Path Algorithm to compute optimal load blocks from a weighted cost-matrix
#   Returns a sequence of optimal indices for ordered dataset

function BellmanFord(graph, src, dest, blks)
    #Initialise
    n = length(graph[1, :])
    V = fill(Inf, n, blks)
    P = fill(Inf, n, blks)

    for k = 1:blks
        V[1, k] = 0
        P[1, k] = 1
    end

    for j = 2:n
        for k = 1:blks
            V[j, k] = graph[1, j]
            P[j, k] = 1
        end

        for i = 1:j-1
            for k = 2:blks
                if V[i, k-1] + graph[i, j] < V[j, k]
                    V[j, k] = V[i, k-1] + graph[i, j]
                    P[j, k] = i
                end
            end
        end
    end

    seq = Array{Int64}(undef, 1, blks + 1)
    seq[blks+1] = n
    seq[blks] = P[n, blks]

    for i = blks-1:-1:1
        prevState = P[seq[i+1], i]
        seq[i] = prevState
    end

    return seq
end

##  Compute block duration of optimal blocks

function blockDuration(blocks, indices)
    for i = 2:length(indices)
        blocks[i-1] = (indices[i] - indices[i-1]) / 2
    end

end

## Write LB TO CSV

function writeLBtoCSV(filename, filelist, indices, avgLoad, blocklengths, T)

    #Print outputs to new file
    # writefile = joinpath(@__DIR__, ("LDC_" * filename * ".csv"))
    writefile = joinpath(dirrr, ("LDC_" * filename * ".csv"))
    file = open(writefile, "w")

    print(file, "LB_Indices")
    for i= 1:length(avgLoad[1,:])
        print(file, ",Mean_Load_",  filelist[i])
    end
    println(file, ",LB_Duration")

    for i = 1:T
        print(file, string(indices[i]))
        for f= 1:length(filelist)
            print(file, "," * string(avgLoad[i,f]))
        end
        println(file, "," * string(blocklengths[i]))
    end

    close(file)
end

## Delete Leap Year Data
#   Leap year data can cause issue with Standardizing that assumes constant n

function deleteLY(data::Array)
    isLY= length(findall(x -> x == "29/2/*"), data[:,4])
    if isLY > 0
        for i= 1:length(isLY)
            data[]
        end
    end
end

## Visualise the Load block approximation of an input LDC
function plotLB(file::String, loadfile::String, title::String)

        #file = file name
        #n = number of data points in original file that LB defined from (so number of days*48)

        #File name for LB data to be read in
        filename= joinpath(@__DIR__,(file*".csv"))
        data= readdlm(filename, ',')
        data= data[2:end, :]

        #File name for LDC data to be read in
        filename= joinpath(@__DIR__,(loadfile*".csv"))
        loaddata= readdlm(filename, ',')
        loaddata= loaddata[2:end, :]
        numHrs= (length(loaddata[:,2])-1)/2

        #Define a rectangle plotting object for the graph
        #args are width, height, x start and y start
        rectangle(w, h, x, y) = Shape(x.+ [0,w,w,0], y .+ [0,0,h,h])

        #Plot curve
        plotobj= plot(legend=false, margin = 5Plots.mm) #,  background_color = :transparent)                                                            #empty plot object
<<<<<<< HEAD
        plot!(rectangle(data[1,3], data[1,2], 0, 0), color= "blue", opacity = 0.7, linewidth = 0.1)
        for i= 2:length(data[:,1])
            plot!(rectangle(data[i,3], data[i,2] ,sum(data[1:i-1,3]) ,0), opacity=0.7, color= "blue", linewidth = 0.1)
        end

        plot!([x for x=0:0.5:numHrs], loaddata[:,2], linewidth = 3, color= "magenta")
=======
        plot!(rectangle(data[1,3], data[1,2], 0, 0), color= "darkgreen", opacity = 0.5, linewidth = 0.15)
        for i= 2:length(data[:,1])
            plot!(rectangle(data[i,3], data[i,2] ,sum(data[1:i-1,3]) ,0), opacity=0.5, color= "darkgreen", linewidth = 0.15)
        end

        plot!([x for x=0:0.5:numHrs], loaddata[:,2], linewidth = 3, color= "royalblue")
>>>>>>> e36df7e0622d6ef050f4a8658cefa290f9136723

        title!(title)
        yaxis!("Load (MW)")
        xaxis!("Hours Utilised")
        display(plotobj)

        #Save pic
        savefig("LDC.png")
 end

##  Data pre-processing to compute mean demands & cost matrix for DP

function process_file(files::Array)

    # Initialise D, costs for first file
    # n = length(readdlm(joinpath(@__DIR__, (files[1] * ".csv")), ',')[2:end,1])           # Assumption: if multiple files, every file has same n
    n = length(readdlm(joinpath(dirrr, (files[1] * ".csv")), ',')[2:end,1])
    D = zeros(n, n)
    costs = fill(Inf, n, n)

    for f= 1:length(files)
        # filename= joinpath(@__DIR__, (files[f] * ".csv"))
        filename= joinpath(dirrr, (files[f] * ".csv"))
        data = DelimitedFiles.readdlm(filename, ',')
        y = Float64.(data[2:end, 2])
        p = Float64.(data[2:end, 3])

        _inner(D, costs, y, p)
    end
    return D, costs, n
end

function _inner(D, costs, y, p)
    n = length(y)
    Threads.@threads for i = 1:n
        dk = y[i]
        @inbounds begin
            for j = (i+1):n
                dk += y[j]
                db = dk / (j - i + 1)
                D[i, j] += db
                z = 0.0
                for k = i:j
                    z += p[k] * abs(y[k] - db)
                end
                if costs[i,j] == Inf
                    costs[i, j] = z
                else
                    costs[i, j] += z
                end
            end
        end
    end
    return
end

##  Main.
#LDC fitting function based on prices
#= Input:
    # files: array of the names of the .csv files to be read in the same directory as this function
    # T= number of segments desired
    Output: Nothing specifically returned.
            .csv output files corresponding to each input file containing the indices that form
            the beginning of optimal load blocks + mean load in defined blocks
    Define:
    n= number of points
=#

## Orignal Function: Suitable for computing individual sets of load blocks for each input dataset
function LDC(files::Array, T::Int64)

    # Read in and initialise data + arrays needed
    print("\nTask In Progress ....")
    for f = 1:length(files)
        #filename = joinpath(@__DIR__, (files[f] * ".csv"))

        #Data pre-processing to compute optimal Load Blocks
        D, costs, n = process_file([files[f]])

        #Calling functions to compute output
        indices = BellmanFord(costs, 1, n, T)
        avgLoad = meanLoad(D, indices)
        blocklengths = Array{Float64}(undef, 1, T)
        blockDuration(blocklengths, indices)
        writeLBtoCSV(files[f], [files[f]], indices, avgLoad, blocklengths, T)
    end
    println("\rTask Completed             \n")
end

## Standardizing Function: Only returns one set of load blocks for all input datasets.
# IMPORTANT: Ensure each dataset input has same number of data points (Eg. autumn -> autumn etc.)
#            Cases such as daylight savings that cause different dataset sizes will break the code.
#  Assumes files_compute is a subset of files_fit
# (Aka the files in compute are also in fit array)

function LDC_std(files_build::Array, files_fit::Array, T::Int64)

        print("\nTask In Progress ....")

        #Data pre-processing to compute optimal Load Blocks
        D, costs, n = process_file(files_build)

        #Calling functions to compute output
        indices = BellmanFord(costs, 1, n, T)
        avgLoad= zeros(Float64, T, length(files_fit))
        for f= 1:length(files_fit)
            D1, = process_file([files_fit[f]])
            avgLoad[:,f]= meanLoad(D1, indices)
        end

        blocklengths = Array{Float64}(undef, 1, T)
        blockDuration(blocklengths, indices)

        # Added the season type to the name of the CSV output file
        writeLBtoCSV("LB_fitted_Data_"*split(files_fit[1], "_")[1], files_fit, indices, avgLoad, blocklengths, T)

        println("\rTask Completed             \n")
end

global dirrr = @__DIR__
# LDC(["netdemandddtest_Summer_2019_0"], 10)

# plotLB("LDC_netdemand_Summer_2019", "netdemand_Summer_2019", "Example Season: 10 Block LDC Approximation")
plotLB("LDC_df_net_00000_Summer", "df_net_00000_Summer", "2019 Summer: 10 Block Stepwise Approximation")
 ##
# Two options:
# First:
# Run script on all CSV files that start with "netdemand" in @__DIR__

#=    files = []
    for file in split.(glob("netdemand*.csv", @__DIR__), "\\")
        push!(files, split(file[end], ".")[1])
    end=#

# Second:
# Specify the names of files to run

#=

global dirrr = "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_important_trials\\Wind_LDC_2019\\Test"
files_build = []
files_fit = []
push!(files_build, ["Summer_Gross"])
push!(files_fit, ["Summer_Gross", "Summer_Net"])

push!(files_build,["Autumn_Gross"])
push!(files_fit, ["Autumn_Gross", "Autumn_Net"])

push!(files_build,["Winter_Gross"])
push!(files_fit, ["Winter_Gross", "Winter_Net"])

push!(files_build,["Spring_Gross"])
push!(files_fit, ["Spring_Gross", "Spring_Net"])

for i in 1:4
    @time begin
    println("FIT:    ", files_build[i])
    println("BUILD:   ", files_fit[i])
    LDC_std(files_build[i], files_fit[i], 10)
    end
end

# dir = []
#
# push!(dir, "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_important_trials\\Solar_blocksizes_2019_Ver2\\Inv_0")
# push!(dir, "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_important_trials\\Solar_blocksizes_2019_Ver2\\Inv_1000")
# push!(dir, "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_important_trials\\Solar_blocksizes_2019_Ver2\\Inv_2000")
# push!(dir, "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_important_trials\\Solar_blocksizes_2019_Ver2\\Inv_4000")
#
# files_build = []
# files_fit = []
#
# push!(files_build, ["netdemanddeductpreSort_Summer_2019"])
# push!(files_fit, ["netdemanddeductpreSort_Summer_2019", "netdemanddeductpostSort_Summer_2019", "netdemandddtest_Summer_2019", "netdemandddtest_Summer_2019_SOLAR"])
#
# push!(files_build,["netdemanddeductpreSort_Autumn_2019"])
# push!(files_fit, ["netdemanddeductpreSort_Autumn_2019", "netdemanddeductpostSort_Autumn_2019", "netdemandddtest_Autumn_2019", "netdemandddtest_Autumn_2019_SOLAR"])
#
# push!(files_build,["netdemanddeductpreSort_Winter_2019"])
# push!(files_fit, ["netdemanddeductpreSort_Winter_2019", "netdemanddeductpostSort_Winter_2019", "netdemandddtest_Winter_2019", "netdemandddtest_Winter_2019_SOLAR"])
#
# push!(files_build,["netdemanddeductpreSort_Spring_2019"])
# push!(files_fit, ["netdemanddeductpreSort_Spring_2019", "netdemanddeductpostSort_Spring_2019", "netdemandddtest_Spring_2019", "netdemandddtest_Spring_2019_SOLAR"])

# files_build = ["netdemanddeductpreSort_Summer_2019"]
# files_fit= ["netdemanddeductpreSort_Summer_2019", "netdemanddeductpostSort_Summer_2019", "netdemandddtest_Summer_2019", "netdemandddtest_Summer_2019_SOLAR"]
#
# files_build = ["netdemanddeductpreSort_Autumn_2019"]
# files_fit= ["netdemanddeductpreSort_Autumn_2019", "netdemanddeductpostSort_Autumn_2019", "netdemandddtest_Autumn_2019", "netdemandddtest_Autumn_2019_SOLAR"]
#
# files_build = ["netdemanddeductpreSort_Winter_2019"]
# files_fit= ["netdemanddeductpreSort_Winter_2019", "netdemanddeductpostSort_Winter_2019", "netdemandddtest_Winter_2019", "netdemandddtest_Winter_2019_SOLAR"]
#
# files_build = ["netdemanddeductpreSort_Spring_2019"]
# files_fit= ["netdemanddeductpreSort_Spring_2019", "netdemanddeductpostSort_Spring_2019", "netdemandddtest_Spring_2019", "netdemandddtest_Spring_2019_SOLAR"]
#
# dir = []
# push!(dir, "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_important_trials\\Solar_standardised_blocks_2019\\Summer")
# push!(dir, "C:\\Users\\calro\\Documents\\Academic\\JuDGE project\\LDC_important_trials\\Solar_standardised_blocks_2019\\Autumn")
#
# files_build = ["netdemanddeductpostSort_Summer_2019_0","netdemanddeductpostSort_Autumn_2019_0"]
# files_fit = []
#
#
#
# for j in 1:2
#     files_fit = []
#     global dirrr = dir[j]
#
#     for file in split.(glob("netdemand*.csv", dirrr), "\\")
#         push!(files_fit, split(file[end], ".")[1])
#     end
#
#     println("FIT:    ", files_build[j])
#     println("BUILD:   ", files_fit)
#     @time begin
#     LDC_std([files_build[j]], files_fit, 10)
#     end
# end

# for j in 1:2
#     global dirrr = dir[j]
#     for i in 1:2
#         @time begin
#         LDC_std(files_build[i], files_fit[i], 10)
#         end
#     end
# end

# for file in split.(glob("netdemand*.csv", dirrr), "\\")
#     push!(files_fit, split(file[end], ".")[1])
# end

# println("Files to be processed:\n", files)
# LDC(files, 10)                            #Change # of load blocks desired
# @time begin
# LDC_std(files_build, files_fit ,10)
# end

#=
TO DO: For LDC_Std, consider how leap years/different season with different n
       prevent code from working (And a roundabout for that)

       Change the LDC_Std to compute only one set of LB
       for all input datasets, then fit the LB widths to
       input data and find mean demand for each one

       => 2 arrays, files to build load blocks on, then fit on all second ones
       => 2 functions, indices passe into the next one
       -> printing into each csv
       #Delete leap years
=#
=#

max_stars_repo_name: calroughan/EPSILON

max_issues_repo_name: calroughan/EPSILON

max_forks_repo_name: calroughan/EPSILON


text:
#!/usr/bin/env julia

# Lists the temporal extents of give RAW or EVR files

using SimradRaw
using EchoviewEvr
using Filetimes
using Base.Iterators

function rawstartend(rawfilename)
    starttime = 0
    endtime = 0
    open(rawfilename, "r") do f
        while !eof(f)
            dgheader= nothing
            body = nothing
            try
                dgheader, body = readencapsulateddatagram(f, datagramreader=readdatagramheader)
            catch ex
                warn("$ex on $f")
                break
            end
            t = SimradRaw.filetime(dgheader.datetime)
            if starttime == 0
                starttime = t
            end
            endtime = t
        end
    end
    return starttime, endtime
end

function evrstartend(evrfilename)
    ps = polygons(evrfilename)

    xs = flatten([x for (x,y) in ps])

    return minimum(xs), maximum(xs)
end

function main(args)

    if length(args) == 0
        x = now()
        println("$x\t$(Filetimes.filetime(x))")
    else
        for arg in args
            if endswith(uppercase(arg), ".RAW")
                starttime, endtime = rawstartend(arg)
                println("$arg\t$starttime\t$endtime\t$(datetime(starttime))\t$(datetime(endtime))")
                continue
            end
            if endswith(uppercase(arg), ".EVR")
                starttime, endtime = evrstartend(arg)
                println("$arg\t$starttime\t$endtime\t$(datetime(starttime))\t$(datetime(endtime))")
                continue
            end

            y = tryparse(Int64,arg)
            if isnull(y)
                println("$arg\t$(Filetimes.filetime(arg))")
            else
                println("$(datetime(get(y)))\t$arg")
            end
        end
    end

end

main(ARGS)

ext: jl


text:
#
# Collection of step sizes
#
export ConstantStepsize, DecreasingStepsize
export Linesearch, ArmijoLinesearch
export getStepsize!, getInitialStepsize, getLastStepsize
#
# Simple ones
#
"""
    ConstantStepsize <: Stepsize

A functor that always returns a fixed step size.

# Fields
* `length` – constant value for the step size.

# Constructor

    ConstantStepSize(s)

initialize the stepsie to a constant `s`
"""
mutable struct ConstantStepsize <: Stepsize
    length::Float64
    ConstantStepsize(s::Real) = new(s)
end
(s::ConstantStepsize)(p::P,o::O,i::Int) where {P <: Problem, O <: Options} = s.length
getInitialStepsize(s::ConstantStepsize) = s.length

@doc doc"""
    DecreasingStepsize()

A functor that represents several decreasing step sizes

# Fields
* `length` – (`1`) the initial step size $l$.
* `factor` – (`1`) a value $f$ to multiply the initial step size with every iteration
* `subtrahend` – (`0`) a value $a$ that is subtracted every iteration
* `exponent` – (`1`) a value $e$ the current iteration numbers $e$th exponential
  is taken of

In total the complete formulae reads for the $i$th iterate as

$ s_i = \frac{(l-i\cdot a)f^i}{i^e}$

and hence the default simplifies to just $ s_i = \frac{l}{i} $

# Constructor

    ConstantStepSize(l,f,a,e)

initialiszes all fields above, where none of them is mandatory.
"""
mutable struct DecreasingStepsize <: Stepsize
    length::Float64
    factor::Float64
    subtrahend::Float64
    exponent::Float64
    DecreasingStepsize(l::Real=1.0, f::Real=1.0, a::Real=0.0, e::Real=1.0) = new(l,f,a,e)
end
(s::DecreasingStepsize)(p::P,o::O,i::Int) where {P <: Problem, O <: Options} = (s.length - i*s.subtrahend)*(s.factor^i)/(i^(s.exponent))
getInitialStepsize(s::DecreasingStepsize) = s.length

#
# Linesearch
#
"""
    Linesearch <: Stepsize

An abstract functor to represent line search type step size deteminations, see
[`Stepsize`](@ref) for details. One example is the [`ArmijoLinesearch`](@ref)
functor.

Compared to simple step sizes, the linesearch functors provide an interface of
the form `(p,o,i,η) -> s` with an additional (but optional) fourth parameter to
proviade a search direction; this should default to something reasonable, e.g.
the negative gradient.
"""
abstract type Linesearch <: Stepsize end

@doc doc"""
    ArmijoLineseach <: Linesearch

A functor representing Armijo line seach including the last runs state, i.e. a
last step size.

# Fields
* `initialStepsize` – (`1.0`) and initial step size
* `retraction` – ([`exp`](@ref Manopt.exp)) the rectraction used in line search
* `contractionFactor` – (`0.95`) exponent for line search reduction
* `sufficientDecrease` – (`0.1`) gain within Armijo's rule
* `lastStepSize` – (`initialstepsize`) the last step size we start the search with

# Constructor

    ArmijoLineSearch()

with the Fields above in their order as optional arguments.

This method returns the functor to perform Armijo line search, where two inter
faces are available:
* based on a tuple `(p,o,i)` of a [`GradientProblem`](@ref) `p`, [`Options`](@ref) `o`
  and a current iterate `i`.
* with `(M,x,F,∇Fx[,η=-∇Fx]) -> s` where [`Manifold`](@ref) `M`, a current
  [`MPoint`](@ref) `x` a function `F`, that maps from the manifold to the reals,
  its gradient (a [`TVector`](@ref)) `∇F`$=\nabla F(x)$ at  `x` and an optional
  search direction [`TVector`](@ref) `η-∇F` are the arguments.
"""
mutable struct ArmijoLinesearch <: Linesearch
    initialStepsize::Float64
    retraction::Function
    contractionFactor::Float64
    sufficientDecrease::Float64
    stepsizeOld::Float64
    ArmijoLinesearch(
        s::Float64=1.0,
        r::Function=exp,
        contractionFactor::Float64=0.95,
        sufficientDecrease::Float64=0.1) = new(s, r, contractionFactor, sufficientDecrease,s)
end
function (a::ArmijoLinesearch)(p::P,o::O,i::Int, η::T=-getGradient(p,o.x)) where {P <: GradientProblem{mT} where mT <: Manifold, O <: Options, T <: TVector}
    a(p.M, o.x, p.costFunction, getGradient(p,o.x), η)
end
function (a::ArmijoLinesearch)(M::mT, x::P, F::Function, ∇F::T, η::T=-∇F) where {mT <: Manifold, P <: MPoint, T <: TVector}
    # for local shortness
    s = a.stepsizeOld
    retr = a.retraction
    f0 = F(x)
    xNew = retr(M,x,s*η)
    fNew = F(xNew)
    while fNew < f0 + a.sufficientDecrease*s*dot(M, x, η, ∇F) # increase
        xNew = retr(M,x,s*η)
        fNew = F(xNew)
        s = s/a.contractionFactor
    end
    s = s*a.contractionFactor # correct last
    while fNew > f0 + a.sufficientDecrease*s*dot(M, x, η, ∇F) # decrease
        s = a.contractionFactor * s
        xNew = retr(M,x,s*η)
        fNew = F(xNew)
    end
    a.stepsizeOld = s
    return s
end
getInitialStepsize(a::ArmijoLinesearch) = a.initialStepsize

#
# Access functions
#
@traitfn getStepsize!(p::P, o::O,vars...) where {P <: Problem, O <: Options; IsOptionsDecorator{O}} = getStepsize!(p, o.options,vars...)
@traitfn getStepsize!(p::P, o::O,vars...) where {P <: Problem, O <: Options; !IsOptionsDecorator{O}} = o.stepsize(p,o,vars...)

@traitfn getInitialStepsize(p::P,o::O,vars...) where {P <: Problem, O <: Options; !IsOptionsDecorator{O}} = getInitialStepsize(o.stepsize)
@traitfn getInitialStepsize(p::P, o::O,vars...) where {P <: Problem, O <: Options; IsOptionsDecorator{O}} = getInitialStepsize(p, o.options)

@traitfn getLastStepsize(p::P, o::O,vars...) where {P <: Problem, O <: Options; IsOptionsDecorator{O}} = getLastStepsize(p, o.options,vars...)
@traitfn getLastStepsize(p::P, o::O,vars...) where {P <: Problem, O <: Options; !IsOptionsDecorator{O}} = getLastStepsize(p,o,o.stepsize,vars...)

getLastStepsize(p::P,o::O,s::S,vars...) where {P <: Problem, O <: Options,S <: Stepsize} = s(p,o,vars...)
getLastStepsize(p::P,o::O,s::ArmijoLinesearch,vars...) where {P <: Problem, O <: Options} = s.stepsizeOld


text:
function transformoutline(outline::ScenarioOutline)
    # The examples are in a multidimensional array.
    # The size of dimension 1 is the number of placeholders.
    # The size of dimension 2 is the number of examples.
    [interpolatexample(outline, outline.examples[:,exampleindex])
     for exampleindex in 1:size(outline.examples, 2)]
end

function interpolatexample(outline::ScenarioOutline, example::Vector{String})
    placeholders_kv = ["<$(outline.placeholders[i])>" => example[i] for i in 1:length(example)]
    placeholders = Dict{String, String}(placeholders_kv...)

    fromplaceholders = x -> placeholders[x]
    steps = [interpolatestep(step, fromplaceholders) for step in outline.steps]

    Scenario(outline.description, outline.tags, steps)
end

interpolatestep(step::Given, fromplaceholders::Function) = Given(interpolatesteptext(step.text, fromplaceholders);
                                                                 block_text=interpolatesteptext(step.block_text, fromplaceholders))
interpolatestep(step::When, fromplaceholders::Function) = When(interpolatesteptext(step.text, fromplaceholders);
                                                               block_text=interpolatesteptext(step.block_text, fromplaceholders))
interpolatestep(step::Then, fromplaceholders::Function) = Then(interpolatesteptext(step.text, fromplaceholders);
                                                               block_text=interpolatesteptext(step.block_text, fromplaceholders))

interpolatesteptext(text::String, fromplaceholders::Function) = replace(text, r"<[^>]*>", fromplaceholders)


text:
# ---
# jupyter:
#   jupytext:
#     formats: jl,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # PSI4 Acquisition

# ### Easy choice: clone and compile
#
# * Clone from [GitHub](https://github.com/psi4/psi4)
# * Compile via [CMake](https://github.com/psi4/psi4/blob/master/CMakeLists.txt#L16-L123)
# * Build FAQ and guide [here](http://psicode.org/psi4manual/master/build_faq.html)

# ### Easier choice: conda
#
# 0. Get [miniconda](https://conda.io/miniconda.html) or [anaconda](https://www.continuum.io/downloads) installer or script
# 0. Run the installer or "bash" the script on the command line, accepting its license and allowing it to add its installation directory to your `PATH` environment variable: _e.g._, ``bash Anaconda3-4.3.0-Linux-x86_64.sh``
# 0. Create a PSI4 environment named "p4env". Until the 1.1 release, instructions may be found [at the end of this thread](https://github.com/psi4/psi4/issues/466#issuecomment-272589229)
# 0. Activate the environment: ``source activate p4env``
# 0. See [guide](http://psicode.org/psi4manual/master/build_planning.html#faq-runordinarymodule) for any configuration trouble

# ### Easiest choice: conda installer
#
# * **Not available until 1.1 release**
# * ~~Get Psi4+miniconda installer script from [psicode](http://psicode.org/downloads.html) and follow its directions~~

# ### Online choice: binder button
# * **Not available until 1.1 release**

# # PSI4 Boilerplate
#
#

# ### Eponymous Julia module imports

using PyCall: pyimport

psi4 = pyimport("psi4")
np   = pyimport("numpy")

# ### Direct output and scratch
#
# * Output goes to file ``output.dat``
# * Boolean directs overwriting (`true`) rather than appending (`false`).
# * Optionally, redirect scratch away from ``/tmp`` to existing, writable directory

# +
psi4.set_output_file("output.dat", true)

# optional
#psi4.core.IOManager.shared_object().set_default_path("/scratch")
# -

# ### Set memory limits
#
# * Give 500 Mb of memory to PSI4
# * Give 2 Gb of memory for NumPy arrays (quantity for Psi4Julia project, *not* passed to NumPy)
# * Sum of these two should nowhere approach the RAM of your computer

psi4.set_memory(Int(5e8))
numpy_memory = 2

# ### Molecule and Basis
#
# * Covered in detail in subsequent tutorials. This is the quick reference
# * Running _without_ symmetry recommended in Psi4NumPy for simplest expressions

# +
psi4.geometry("""
O 0.0 0.0 0.0
H 1.0 0.0 0.0
H 0.0 1.0 0.0
symmetry c1
""")

psi4.set_options(Dict("basis" => "cc-pvdz"))


text:
function combine_decomposition!(db::SQLite.DB, id::Int, name::AbstractString, scenario::Union{Int, AbstractString},
                                clique_path::AbstractString, cliquetree_path::AbstractString, nb_added_edge_dec::Int,
                                id_1::Int, graph_path_1::AbstractString,
                                id_2::Int, graph_path_2::AbstractString;
                                how::AbstractString, extension_alg::AbstractString, rng)
    println("Combining decompositions ($name, $scenario) (in1=$id_1, in2=$id_2)")

    # Retrieve the graphs
    g1 = loadgraph(graph_path_1)
    g2 = loadgraph(graph_path_2)

    # Combine
    g3 = combine_graph(g1, g2; how=how, extension_alg=extension_alg)

    # Extract features
    clique = maximal_cliques(g3)
    cliquetree, nb_lc = get_cliquetree_array(clique)
    features = get_features_graph(g3)
    features["nb_lc"] = nb_lc
    features["nb_added_edge_dec"] = nb_added_edge_dec + (ne(g3) - ne(g1))
    merge!(features, get_clique_features(clique))

    # Save cliques, cliquetree and graph
    cliques_path = dirname(clique_path)
    cliquetrees_path = dirname(cliquetree_path)
    graphs_path = dirname(graph_path_1)
    uuid = uuid1(rng)
    clique_path = joinpath(cliques_path, "$(name)_$(scenario)_$(uuid)_clique.csv")
    cliquetree_path = joinpath(cliquetrees_path, "$(name)_$(scenario)_$(uuid)_cliquetree.csv")
    graph_path = joinpath(graphs_path, "$(name)_$(scenario)_$(uuid)_graph.lgz")
    save_clique(clique, clique_path)
    save_cliquetree(cliquetree, cliquetree_path)
    savegraph(graph_path, g3)

    # Other columns
    date = Dates.now()

    features = Dict(Symbol(k) => v for (k, v) in features)
    insert_decomposition!(db, id, uuid, name, scenario, "combine", extension_alg, date, clique_path, cliquetree_path, graph_path; features...)

    # Insert in merge table
    out_id = DBInterface.execute(db, "SELECT id FROM decompositions WHERE uuid = '$uuid'") |> DataFrame
    out_id = out_id[1, :id]
    insert_combination!(db, id_1, id_2, out_id, how, extension_alg)
end

function combine_decomposition_dfrow!(db::SQLite.DB, row; how, extension_alg, rng)
    combine_decomposition!(db, row["origin_id"], row["origin_name"], row["origin_scenario"],
                           row["clique_path"], row["cliquetree_path"], row["nb_added_edge_dec"],
                           row["id"], row["graph_path"],
                           row["id_1"], row["graph_path_1"];
                           how=how, extension_alg=extension_alg, rng=rng)
end

function combine_decompositions!(db::SQLite.DB; how::AbstractString, extension_alg::AbstractString,
                                 min_nv=typemin(Int), max_nv=typemax(Int), subset=nothing, exclude=["combine"], rng=MersenneTwister(42))
    query = """
    SELECT d1.origin_id, d1.id, d1.origin_name, d1.origin_scenario, d1.graph_path, d1.clique_path, d1.cliquetree_path, d1.nb_added_edge_dec, d1.extension_alg,
           d2.id, d2.origin_name, d2.origin_scenario, d2.graph_path, d2.extension_alg
    FROM decompositions AS d1 CROSS JOIN decompositions AS d2
    WHERE d1.origin_name = d2.origin_name AND d1.origin_scenario = d2.origin_scenario AND d1.id > d2.id
    """
    if !isnothing(subset)
        query *= " AND d1.origin_id IN ($(join(subset, ',')))"
    end
    for e in exclude
        query *= " AND d1.extension_alg <> '$e' AND d2.extension_alg <> '$e'"
    end
    results = DBInterface.execute(db, query) |> DataFrame
    println("Combine config: how=$how, extension_alg=$extension_alg, min_nv=$min_nv, max_nv=$max_nv, exlucde=$exclude, rng=$rng")
    println("subset\n$subset\nend subset")
    combine_function!(row) = combine_decomposition_dfrow!(db, row; how=how, extension_alg=extension_alg, rng=rng)
    combine_function!.(eachrow(results[!, ["origin_id", "origin_name", "origin_scenario", "clique_path", "cliquetree_path", "nb_added_edge_dec",
                                           "id", "graph_path", "id_1", "graph_path_1"]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
#################### Model Graph ####################

function ModelGraph(m::Model)
  allkeys = keys(m, :all)
  g = DiGraph(length(allkeys))
  lookup = Dict(allkeys[i] => i for i in 1:length(allkeys))
  for key in keys(m)
    node = m[key]
    if isa(node, AbstractDependent)
      for src in node.sources
        add_edge!(g, lookup[src], lookup[key])
      end
    end
  end
  ModelGraph(g, allkeys)
end

#################### Display ####################

function draw(m::Model; filename::AbstractString="")
  dot = graph2dot(m)
  if length(filename) == 0
    print(dot)
  else
    if something(findfirst(isequal('.'), filename), 0) == 0
      filename = string(filename, ".dot")
    end
    f = open(filename, "w")
    write(f, dot)
    close(f)
  end
end

graph(m::Model) = ModelGraph(m)

function graph2dot(m::Model)
  dag = ModelGraph(m)
  io = IOBuffer()
  write(io, "digraph MambaModel {\n")
  deps = keys(m, :dependent)
  for v in vertices(dag.graph)
    attr = Tuple{AbstractString, AbstractString}[]
    vkey = dag.keys[v]
    if vkey in deps
      node = m[vkey]
      if isa(node, AbstractLogical)
        push!(attr, ("shape", "diamond"))
      elseif isa(node, AbstractStochastic)
        push!(attr, ("shape", "ellipse"))
      end
      if isempty(node.monitor)
        push!(attr, ("style", "filled"),
                    ("fillcolor", "gray85"))
      end
    else
      push!(attr, ("shape", "box"),
                  ("style", "filled"),
                  ("fillcolor", "gray85"))
    end
    write(io, "\t\"")
    write(io, vkey)
    write(io, "\" [")
    write(io, join(map(x -> "$(x[1])=\"$(x[2])\"", attr), ", "))
    write(io, "];\n")
    for t in outneighbors(dag.graph, v)
      write(io, "\t\t\"")
      write(io, vkey)
      write(io, "\" -> \"")
      write(io, dag.keys[t])
      write(io, "\";\n")
     end
  end
  write(io, "}\n")
  String(take!(copy(io)))
end

#################### Auxiliary Functions ####################

function any_stochastic(dag::ModelGraph, v::Int, m::Model)
  found = false
  for t in outneighbors(dag.graph, v)
    tkey = dag.keys[t]
    if isa(m[tkey], AbstractStochastic) || any_stochastic(dag, t, m)
      found = true
      break
    end
  end
  found
end

function gettargets(dag::ModelGraph, v::Int, terminalkeys::Vector{Symbol})
  values = Symbol[]
  for t in outneighbors(dag.graph, v)
    tkey = dag.keys[t]
    push!(values, tkey)
    if !(tkey in terminalkeys)
      values = union(values, gettargets(dag, t, terminalkeys))
    end
  end
  values
end

function tsort(m::Model)
  dag = ModelGraph(m)
  dag.keys[topological_sort_by_dfs(dag.graph)]
end


text:
# include -
include("Include.jl")

# setup the calculation -
number_of_samples = 5
number_of_parameters = 153

# load up the parameters sets -
ensemble_array = zeros(number_of_parameters,1);
for sample_index = 1:number_of_samples

  filename = "./parameter_best_v1.dat."*string(sample_index)
  local_array = readdlm(filename)

  ensemble_array = hcat(ensemble_array,local_array);
end

# cut leading col -
ensemble_array = ensemble_array[:,2:end]

# Last ... run the model w/these sets, save the data -
for parameter_set_index = 1:number_of_samples

  # sample the ensemble_array -
  parameter_array = ensemble_array[:,parameter_set_index]

  # Load the data dictionary -
  data_dictionary = DataDictionary(0.0,0.0,0.0);

  # Update data dictionary to match new parameters before calculating obj
  parameter_mapping_array = data_dictionary["parameter_name_mapping_array"]
  for index = 1:length(parameter_mapping_array)
    if index <= data_dictionary["number_of_binding"]
      data_dictionary["binding_parameter_dictionary"][parameter_mapping_array[index]] = parameter_array[index]
    else
      data_dictionary["control_parameter_dictionary"][parameter_mapping_array[index]] = parameter_array[index]
    end
  end

  # Setup the time -
  time_start = 0.0
  time_stop = 120.0
  time_step_size = 0.01

  # Run the simulation -
  (T,X) = Simulation(time_start,time_stop,time_step_size,data_dictionary);

  # dump data to disk -
  local_data = [T X];
  data_filename = "./best_ensemble/sim_data_PSI_"*string(parameter_set_index)*".dat"
  writedlm(data_filename,local_data);
end


text:
import Base.sqrt

function power_of_2_gt(n::Int)
    Int(2^ceil(log(2,n)))
end

struct ScaledGeoMeanAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr, AbstractExpr}
    size::Tuple{Int, Int}

    function ScaledGeoMeanAtom(x::AbstractExpr)
        children = (x)
        return new(:scaledgeomean, hash(children), children, x.size)
    end
end

function sign(q::ScaledGeoMeanAtom)
    return Positive()
end

function monotonicity(q::ScaledGeoMeanAtom)
    return (Nondecreasing(),)
end

function curvature(q::ScaledGeoMeanAtom)
    return ConcaveVexity()
end

function evaluate(q::ScaledGeoMeanAtom)
    nbar = power_of_2_gt(length(q.children))
    return prod(evaluate(q.children[1]))^(1/nbar)
end

function conic_form!(q::ScaledGeoMeanAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, q)
        t = Variable()
        qol_objective = conic_form!(t, unique_conic_forms)
        x, y = q.children
        conic_form!(SOCElemConstraint(y + x, y - x, 2 * t), unique_conic_forms)
        conic_form!(x >= 0, unique_conic_forms)
        conic_form!(y >= 0, unique_conic_forms)
        cache_conic_form!(unique_conic_forms, q, qol_objective)
    end
    return get_conic_form(unique_conic_forms, q)
end

function scaledgeomean(x::AbstractExpr)
    if length(x) > 2
        return ScaledGeoMeanAtom(x)
    elseif length(x) == 2
        return geomean(x[1],x[2])
    else
        return x
    end
end


text:
export
    PDB,
    PDBParseError,
    downloadpdb,
    spaceatomname,
    pdbline,
    writepdb

"Protein Data Bank (PDB) file format."
immutable PDB <: Bio.IO.FileFormat end

"Error arising from parsing a Protein Data Bank (PDB) file."
type PDBParseError <: Exception
    message::String
    line_number::Int
    line::String
end

function Base.showerror(io::IO, e::PDBParseError)
    return print(io,
            e.message,
            " at line ",
            e.line_number,
            " of file:\n",
            e.line)
end

"""
Download a Protein Data Bank (PDB) file or biological assembly from the RCSB
PDB. By default downloads the PDB file; if `ba_number` is set the biological
assembly with that number will be downloaded.
"""
function downloadpdb(pdbid::AbstractString,
                    out_filepath::AbstractString="$pdbid.pdb";
                    ba_number::Integer=0)
    # Check PDB ID is 4 characters long and only consits of alphanumeric characters
    if length(pdbid) != 4 || ismatch(r"[^a-zA-Z0-9]", pdbid)
        throw(ArgumentError("Not a valid PDB ID: \"$pdbid\""))
    end
    if ba_number == 0
        download("http://www.rcsb.org/pdb/files/$pdbid.pdb", out_filepath)
    else
        # Will download error page if ba_number is too high
        download("http://www.rcsb.org/pdb/files/$pdbid.pdb$ba_number", out_filepath)
    end
end

function Base.read(input::IO,
            ::Type{PDB};
            structure_name::AbstractString="",
            remove_disorder::Bool=false,
            read_std_atoms::Bool=true,
            read_het_atoms::Bool=true)
    # Define ProteinStructure and add to it incrementally
    struc = ProteinStructure(structure_name)
    struc[1] = Model(1, struc)
    # Entries outside of a MODEL/ENDMDL block are added to model 1
    curr_model = 1
    line_n = 0
    for line in eachline(input)
        line_n += 1
        # Read ATOM and HETATM records as required
        if (read_std_atoms && startswith(line, "ATOM  ")) ||
                (read_het_atoms && startswith(line, "HETATM"))
            unsafe_addatomtomodel!(
                    struc[curr_model],
                    AtomRecord(line, line_n), remove_disorder=remove_disorder)
        # Read MODEL record
        elseif startswith(line, "MODEL ")
            try
                curr_model = parse(Int, line[11:14])
            catch
                throw(PDBParseError(
                    "Could not read model serial number", line_n, line))
            end
            # Create model if required
            if !haskey(models(struc), curr_model)
                struc[curr_model] = Model(curr_model, struc)
            end
        # Read ENDMDL record
        elseif startswith(line, "ENDMDL")
            curr_model = 1
        end
    end
    # Remove any models that were not added to
    for model in struc
        if countchains(model) == 0
            delete!(models(struc), modelnumber(model))
        end
    end
    # Generate lists for iteration
    fixlists!(struc)
    return struc
end

function Base.read(filepath::AbstractString,
            ::Type{PDB};
            structure_name::AbstractString=splitdir(filepath)[2],
            kwargs...)
    open(filepath, "r") do input
        read(input, PDB; structure_name=structure_name, kwargs...)
    end
end

# Constructor from PDB ATOM/HETATM line
AtomRecord(pdb_line::String, line_n::Integer=1) = AtomRecord(
    pdb_line[1] == 'H', # This assumes the line has already been checked as an ATOM/HETATM record
    parsestrict(pdb_line, 7, 11, Int, "Could not read atom serial number", line_n),
    parsestrict(pdb_line, 13, 16, String, "Could not read atom name", line_n),
    parsestrict(pdb_line, 17, 17, Char, "Could not read alt loc identifier", line_n),
    parsestrict(pdb_line, 18, 20, String, "Could not read residue name", line_n),
    parsestrict(pdb_line, 22, 22, Char, "Could not read chain ID", line_n),
    parsestrict(pdb_line, 23, 26, Int, "Could not read residue number", line_n),
    parsestrict(pdb_line, 27, 27, Char, "Could not read insertion code", line_n),
    [
        parsestrict(pdb_line, 31, 38, Float64, "Could not read x coordinate", line_n),
        parsestrict(pdb_line, 39, 46, Float64, "Could not read y coordinate", line_n),
        parsestrict(pdb_line, 47, 54, Float64, "Could not read z coordinate", line_n)
    ],
    parselenient(pdb_line, 55, 60, Float64, 1.0),
    parselenient(pdb_line, 61, 66, Float64, 0.0),
    parselenient(pdb_line, 77, 78, String, "  "),
    parselenient(pdb_line, 79, 80, String, "  ")
)

"Parse columns from a line and return the value or throw a `PDBParseError`."
function parsestrict(line::String,
                    col_1::Integer,
                    col_2::Integer,
                    out_type::Type,
                    error_message::AbstractString,
                    line_n::Integer)
    try
        return parsevalue(line, col_1, col_2, out_type)
    catch
        throw(PDBParseError(error_message, line_n, line))
    end
end

"Parse columns from a line and return the value or a default value."
function parselenient(line::String,
                    col_1::Integer,
                    col_2::Integer,
                    out_type::Type,
                    default)
    try
        return parsevalue(line, col_1, col_2, out_type)
    catch
        return default
    end
end

"Parse columns from a line into a given type."
function parsevalue(line::String,
                    col_1::Integer,
                    col_2::Integer,
                    out_type::Type)
    try
        if out_type == Int
            return parse(Int, line[col_1:col_2])
        elseif out_type == Float64
            return parse(Float64, line[col_1:col_2])
        elseif out_type == String
            return line[col_1:col_2]
        elseif out_type == Char
            return line[col_1]
        else
            error()
        end
    catch
        error("Could not parse to desired type")
    end
end

"""
Form a string of a certain length from a value by adding spaces to the left.
Throws an error if the value is too long.
"""
function spacestring(val_in, new_length::Integer)
    string_out = string(val_in)
    if length(string_out) > new_length
        throw(ArgumentError("Cannot fit value \"$string_out\" into $new_length space(s)"))
    end
    return lpad(string_out, new_length)
end

"""
Space an `Atom` name such that the last element letter (generally) appears in
the second column. If the `element` property of the `Atom` is set it is used to
get the element, otherwise the name starts from the second column where
possible. This function is generally not required as spacing is recorded when
atom names are read in from a Protein Data Bank (PDB) file.
"""
function spaceatomname(at::Atom)
    at_name = atomname(at, strip=false)
    chars = length(at_name)
    if chars == 4
        return at_name
    end
    strip_el = element(at, strip=true)
    if chars > 4
        throw(ArgumentError("Atom name is greater than four characters: \"$at_name\""))
    end
    # In the absence of the element, the first index goes in column two
    if strip_el == "" || findfirst(at_name, strip_el[1]) == 0
        cent_ind = 1
    # The last letter of the element goes in column two where possible
    else
        cent_ind = findfirst(at_name, strip_el[1]) + length(strip_el) - 1
    end
    if cent_ind > 2
        throw(ArgumentError("Atom name is too long to space correctly: \"$at_name\""))
    end
    if cent_ind == 1 && chars < 4
        out_string = " $at_name"
    else
        out_string = "$at_name"
    end
    return rpad(out_string, 4)
end

"""
Form a Protein Data Bank (PDB) format ATOM or HETATM record from an `Atom` or
`AtomRecord`.
"""
function pdbline(at::Atom)
    return (ishetero(at) ? "HETATM" : "ATOM  ") *
            spacestring(serial(at), 5) *
            " " *
            spaceatomname(at) *
            string(altlocid(at)) *
            spacestring(resname(at), 3) *
            " " *
            string(chainid(at)) *
            spacestring(resnumber(at), 4) *
            string(inscode(at)) *
            "   " *
            # This will throw an error for large coordinate values, e.g. -1000.123
            spacestring(round(x(at), 3), 8) *
            spacestring(round(y(at), 3), 8) *
            spacestring(round(z(at), 3), 8) *
            spacestring(round(occupancy(at), 2), 6) *
            # This will throw an error for large temp facs, e.g. 1000.12
            spacestring(round(tempfactor(at), 2), 6) *
            "          " *
            spacestring(element(at), 2) *
            spacestring(charge(at), 2)
end

function pdbline(at_rec::AtomRecord)
    return (at_rec.het_atom ? "HETATM" : "ATOM  ") *
            spacestring(at_rec.serial, 5) *
            " " *
            spacestring(at_rec.atom_name, 4) *
            string(at_rec.alt_loc_id) *
            spacestring(at_rec.res_name, 3) *
            " " *
            string(at_rec.chain_id) *
            spacestring(at_rec.res_number, 4) *
            string(at_rec.ins_code) *
            "   " *
            # This will throw an error for large coordinate values, e.g. -1000.123
            spacestring(round(at_rec.coords[1], 3), 8) *
            spacestring(round(at_rec.coords[2], 3), 8) *
            spacestring(round(at_rec.coords[3], 3), 8) *
            spacestring(round(at_rec.occupancy, 2), 6) *
            # This will throw an error for large temp facs, e.g. 1000.12
            spacestring(round(at_rec.temp_factor, 2), 6) *
            "          " *
            spacestring(at_rec.element, 2) *
            spacestring(at_rec.charge, 2)
end

"""
Write a `StructuralElementOrList` to a Protein Data Bank (PDB) format file. Only
ATOM, HETATM, MODEL and ENDMDL records are written - there is no header and
there are no TER records.
Additional arguments are atom selector functions - only atoms that return
`true` from the functions are retained.
"""
function writepdb(output::IO,
                el::Union{ProteinStructure, Vector{Model}},
                atom_selectors::Function...)
    # If there are multiple models, write out MODEL/ENDMDL lines
    if length(el) > 1
        for mod in sort(collect(el))
            println(output, "MODEL     ", spacestring(modelnumber(mod), 4), repeat(" ", 66))
            writepdb(output, mod, atom_selectors...)
            println(output, "ENDMDL$(repeat(" ", 74))")
        end
    # If there is only one model, do not write out MODEL/ENDMDL lines
    elseif isa(el, ProteinStructure)
        writepdb(output, defaultmodel(el), atom_selectors...)
    else
        writepdb(output, el[1], atom_selectors...)
    end
end

function writepdb(output::IO,
                el::Union{Model, Chain, AbstractResidue, Vector{Chain},
                    Vector{AbstractResidue}, Vector{Residue}, Vector{DisorderedResidue}},
                atom_selectors::Function...)
    # Collect residues then expand out disordered residues
    for res in collectresidues(el)
        if isa(res, Residue)
            for at in collectatoms(res, atom_selectors...)
                writepdb(output, at)
            end
        else
            for res_name in resnames(res)
                for at in collectatoms(disorderedres(res, res_name), atom_selectors...)
                    writepdb(output, at)
                end
            end
        end
    end
end

function writepdb(output::IO, at::AbstractAtom, atom_selectors::Function...)
    for atom_record in at
        println(output, pdbline(atom_record))
    end
end

function writepdb{T <: AbstractAtom}(output::IO,
                    ats::Vector{T},
                    atom_selectors::Function...)
    for at in collectatoms(ats, atom_selectors...)
        writepdb(output, at)
    end
end

function writepdb(filepath::AbstractString,
                el::StructuralElementOrList,
                atom_selectors::Function...)
    open(filepath, "w") do output
        writepdb(output, el, atom_selectors...)
    end
end


text:
_fieldtypes(::Type{<:StructArrays.StructArray{<:Any,<:Any,C}}) where {C<:Tuple} = C
_fieldtypes(
    ::Type{<:StructArrays.StructArray{<:Any,<:Any,<:NamedTuple{<:Any,C}}},
) where {C<:Tuple} = C

isimmutableindex(::Type{T}) where {T<:StructArrays.StructArray} =
    _tall(isimmutableindex, _fieldtypes(T))
isimmutablevalue(::Type{T}) where {T<:StructArrays.StructArray} =
    _tall(isimmutablevalue, _fieldtypes(T))
isfrozenindex(::Type{T}) where {T<:StructArrays.StructArray} =
    _tall(isfrozenindex, _fieldtypes(T))
isfrozenvalue(::Type{T}) where {T<:StructArrays.StructArray} =
    _tall(isfrozenvalue, _fieldtypes(T))

ownmutableindex(x::StructArrays.StructArray{T,N}) where {T,N} =
    StructArrays.StructArray{T,N}(_map(ownmutableindex, StructArrays.fieldarrays(x)))


text:
"""Function to import d18O data and process yearmarkers and calculation windows.

Takes the name of a file that is formatted according to the standard format
and converts it to an object to be used later in the model. In doing so, the
function also reads the user-provided yearmarkers in the file and uses them
as a basis for the length of windows used throughout the model. This ensures
that windows are not too short and by default contain at least one year of
growth for modeling.

"""
function data_import(filename::String)
    # Load CSV file
    dat = CSV.read(filename, DataFrame)

    # Replace species column of dummy dataset with dummy yearmarkaers
    dat[:, 3], dat[:, 5] = dat[:, 5], dat[:, 3]
    yearmarkers = rand([true, false], nrow(dat))

    # Check data format
    @assert(ncol(dat) ∈ [3, 5],
        "Wrong number of columns. Abort.")
    if ncol(dat) == 3
        dist = measurement.(dat[:, 1], missing)
        d18o = measurement.(dat[:, 3], missing)
        # yearmarkers = dat[:, 2]
    elseif ncol(dat) == 5
        dist = measurement.(dat[:, 1], dat[:, 2])
        d18o = measurement.(dat[:, 4], dat[:, 5])
        # yearmarkers = dat[:, 3]
    end
    print(yearmarkers[1:5])

    # Check if yearmarkers are correct
    @assert(eltype(yearmarkers) == Bool,
        "Yearmarkers not binary.")
    @assert(count(yearmarker) <= 1,
        "Less than a full year sampled. Abort.")

    modeldata = (; dist, d18o, yearmarkers)

    return modeldata
end


text:
################################################################################

using HiddenMarkovModelReaders
using Test

################################################################################

tests = [
  "distance.jl",
  "setup.jl",
  "process.jl",
]

################################################################################

@testset "HiddenMarkovModelReaders.jl" begin

  for tst ∈ tests
    include(tst)
  end

end

################################################################################

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: appintegrations
using AWS.Compat
using AWS.UUIDs

"""
    CreateEventIntegration()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Creates an EventIntegration, given a specified name, description, and a reference to an Amazon Eventbridge bus in your account and a partner event source that will push events to that bus. No objects are created in the your account, only metadata that is persisted on the EventIntegration control plane.

# Required Parameters
- `EventBridgeBus`: The Eventbridge bus.
- `EventFilter`: The event filter.
- `Name`: The name of the event integration.

# Optional Parameters
- `ClientToken`: A unique, case-sensitive identifier that you provide to ensure the idempotency of the request.
- `Description`: The description of the event integration.
- `Tags`: One or more tags.
"""
create_event_integration(EventBridgeBus, EventFilter, Name; aws_config::AWSConfig=global_aws_config()) = appintegrations("POST", "/eventIntegrations", Dict{String, Any}("EventBridgeBus"=>EventBridgeBus, "EventFilter"=>EventFilter, "Name"=>Name, "ClientToken"=>string(uuid4())); aws_config=aws_config)
create_event_integration(EventBridgeBus, EventFilter, Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("POST", "/eventIntegrations", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EventBridgeBus"=>EventBridgeBus, "EventFilter"=>EventFilter, "Name"=>Name, "ClientToken"=>string(uuid4())), args)); aws_config=aws_config)

"""
    DeleteEventIntegration()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Deletes the specified existing event integration. If the event integration is associated with clients, the request is rejected.

# Required Parameters
- `Name`: The name of the event integration.

"""
delete_event_integration(Name; aws_config::AWSConfig=global_aws_config()) = appintegrations("DELETE", "/eventIntegrations/$(Name)"; aws_config=aws_config)
delete_event_integration(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("DELETE", "/eventIntegrations/$(Name)", args; aws_config=aws_config)

"""
    GetEventIntegration()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Return information about the event integration.

# Required Parameters
- `Name`: The name of the event integration.

"""
get_event_integration(Name; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/eventIntegrations/$(Name)"; aws_config=aws_config)
get_event_integration(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/eventIntegrations/$(Name)", args; aws_config=aws_config)

"""
    ListEventIntegrationAssociations()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Returns a paginated list of event integration associations in the account.

# Required Parameters
- `Name`: The name of the event integration.

# Optional Parameters
- `maxResults`: The maximum number of results to return per page.
- `nextToken`: The token for the next set of results. Use the value returned in the previous response in the next request to retrieve the next set of results.
"""
list_event_integration_associations(Name; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/eventIntegrations/$(Name)/associations"; aws_config=aws_config)
list_event_integration_associations(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/eventIntegrations/$(Name)/associations", args; aws_config=aws_config)

"""
    ListEventIntegrations()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Returns a paginated list of event integrations in the account.

# Optional Parameters
- `maxResults`: The maximum number of results to return per page.
- `nextToken`: The token for the next set of results. Use the value returned in the previous response in the next request to retrieve the next set of results.
"""
list_event_integrations(; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/eventIntegrations"; aws_config=aws_config)
list_event_integrations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/eventIntegrations", args; aws_config=aws_config)

"""
    ListTagsForResource()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Lists the tags for the specified resource.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the resource.

"""
list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    TagResource()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Adds the specified tags to the specified resource.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the resource.
- `tags`: One or more tags.

"""
tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = appintegrations("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Removes the specified tags from the specified resource.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) of the resource.
- `tagKeys`: The tag keys.

"""
untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = appintegrations("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateEventIntegration()

The Amazon AppIntegrations APIs are in preview release and are subject to change. Updates the description of an event integration.

# Required Parameters
- `Name`: The name of the event integration.

# Optional Parameters
- `Description`: The description of the event inegration.
"""
update_event_integration(Name; aws_config::AWSConfig=global_aws_config()) = appintegrations("PATCH", "/eventIntegrations/$(Name)"; aws_config=aws_config)
update_event_integration(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = appintegrations("PATCH", "/eventIntegrations/$(Name)", args; aws_config=aws_config)


text:
"""
    Molecules.function parse_file(file::String; unit=:angstrom, F=Float64, I=Int16)

Reads a xyz file and return a vector o `Atom` objects. Units can be indicated through the keyword argument `unit`.
Coordinates and mass are stored as `F` (`Float64` by default) and atomic number is taken as `I` (`Int16` by default),
where `F` and `I` are keyword arguments.

See also: parse_string
"""
function parse_file(file::String; unit=:angstrom, F=Float64, I=Int16)
    molstring = read(file, String)
    parse_string(molstring::String; unit=unit, F=F, I=I)
end

"""
    Molecules.function parse_string(molstring::String; unit=:angstrom, F=Float64, I=Int16)

Reads a String representing a XYZ file and return a vector o `Atom` objects. Units can be indicated through the keyword argument `unit`.
Coordinates and mass are stored as `F` (`Float64` by default) and atomic number is taken as `I` (`Int16` by default),
where `F` and `I` are keyword arguments.

See also: parse_file
"""
function parse_string(molstring::String; unit=:angstrom, F=Float64, I=Int16)

    # Get a list of Atom objects from String
    if unit == :bohr
        conv = bohr_to_angstrom
    elseif unit == :angstrom
        conv = 1.0
    else
        throw(ArgumentError("unsupported unit requested for XYZ parsing: $unit"))
    end

    # Regex for parsing
    # Match a floating point, e.g. XYZ coordinate or mass, with at least one blank space on the left of it
    re_float = r"\s+([+-]?(?:\d+[.]?\d*|\d*[.]?\d+))"
    # Match an atom id (String for atomic symbol or integer for atomic number) with one or more blank spaces on the left of it
    re_id = r"\s*(\d+|\w{1,2})"

    atoms = Atom{I,F}[]

    # Possible line formats
    # ID Mass X Y Z
    r1 = re_id * re_float * re_float * re_float * re_float
    # ID X Y Z
    r2 = re_id * re_float * re_float * re_float

    line_num = 1
    for line in split(strip(molstring), "\n")

        if contains(line, r1)
            m = match(r1, line)
            # id is taken as a Symbol if a String is given ("H" → :H), otherwise convert to number ("6" → 6)
            id = occursin(r"\d+", m.captures[1]) ? parse(Int64, m.captures[1]) : Symbol(m.captures[1])

            # Try to parse Mass string into a Float
            try
                mass = parse(F, m.captures[2])
            catch ArgumentError
                throw(ArgumentError("Failed to process data in line $line_num:\n $(line)"))
            end
            str_xyz = m.captures[3:5]

        elseif contains(line, r2)
            m = match(r2, line)
            # id is taken as a Symbol if a String is given ("H" → :H), otherwise convert to number ("6" → 6)
            id = occursin(r"\d+", m.captures[1]) ? parse(Int64, m.captures[1]) : Symbol(m.captures[1])

            # Get mass from PeriodicTable
            mass = convert(F, elements[id].atomic_mass / 1u"u")
            str_xyz = m.captures[2:4]
        else
            throw(ArgumentError("Failed to process data in line $line_num:\n $(line)"))
        end

        Z = I(elements[id].number)

        # Convert String vector to Float vector
        xyz = zeros(F, 3)
        try
            xyz .= parse.(F, str_xyz).*conv
        catch ArgumentError
            throw(ArgumentError("Failed to process XYZ coordinates in line $line_num:\n $(m[2:4])"))
        end

        push!(atoms, Atom(Z, mass, xyz))
        line_num += 1
    end

    return atoms
end

"""
    get_xyz(M::Molecule)

Returns a XYZ string in angstrom for the given Molecule.
"""
function get_xyz(M::Molecule)
    molstring = ""
    for A in M
        molstring *= format("{}   {: 15.12f}   {: 15.12f}   {: 15.12f}\n", Molecules.symbol(A), A.xyz...)
    end
    return molstring
end


text:
using BenchmarkTools

A = rand(1000,1000)
B = rand(1000,1000)
C = rand(1000,1000)

function fused(A,B,C)
    A .+ B .+ C
end
@btime fused(A,B,C)
δ = rand(10,10)
D = zeros(1000,1000)
function fused_1!(D,A,B,C)
    D = A .+ B .+ C # equiv to tmp = A .+ B .+ C; D = tmp
end
δ = fused_1!(D,A,B,C)
@btime fused_1!(D,A,B,C)

D = zeros(1000,1000)
function fused_2!(D,A,B,C)
    for j in 1:2000, i in 1:2000
        D[i,j] = A[i,j] + B[i,j] + C[i,j]
    end
    return D
end
δ = fused_2!(D,A,B,C)

D = zeros(1000,1000)
function fused_3!(D,A,B,C)
D .= A .+ B .+ C
end
fused_3!(D,A,B,C)
δ = fused_3!(D,A,B,C)
@btime fused_2!(D,A,B,C)
@btime fused_3!(D,A,B,C)

a = [1,2,3]
function f!(a)
a[1] = 10
end
f!(a)
a # a = [10,2,3]
function g!(a)
a = [1,2]
end
g!(a)
a # a = [10,2,3]

@btime f!(a) # a[1] = 10
@btime g!(a) # a = [1,2]

function g!(a)
    deleteat!(a, 3)
    a[:] = [1,2]
end
g!(a)
a # a = [10,2,3]

function f_alloc(A)
    B = A[1:5,1:5]
end

function f_point(A)
    B = @view A[1:5,1:5]
end

@btime f_alloc(A) # allocate array, 80 ns
@btime f_point(A)

A = rand(1000,1000) # Bind A to 1000x1000 rand matrix
B = f_alloc(A) # Copy 5x5 slice of A to new allocation, B
A .= zeros(1000,1000) # change contents of A
B # B is unchanged

A = rand(1000,1000) # Bind A to new 1000x1000 rand matrix
B = f_point(A) # Bind B to a pointer, which points at 5x5 slice of A
A .= zeros(1000,1000) # change contents of A

function f!(B)
    B[1,2]=10
end

f!(B)

B[1,2]

A[1,2]

N = 1000000
function f(x,y,N)
    z=0
    for i=1:N
        z += x*y
    end
 end

 @btime f(0.0,3.0)

 function f(x,y)
     z=0
     for i=1:10
         z += x*y
     end
  end

  @btime f(2,3)

@code_llvm f(2,5)

@code_llvm f(2.0,5.0)

@code_typed f(2,5.0)
@code_llvm f(2,5.0)

ff(x::Int,y::Int) = 2x + y
ff(x::Float64,y::Float64) = x/y
@show ff(2,5)
@show ff(2.0,5.0)

ff(x::Number,y::Number) = 2x
ff(x::Union{Int,Float64},y::Union{Int,Float64}) = 10 + x

ff(2.0,5)
typeof(2.0) <: Int

ff(x::Float64,y::Number) = 5x + 2y
ff(x::Number,y::Int) = x - y
ff(2.0,5)

function h()
    a = [1,2,3]
end
@btime h()

function hh()
    a = [0.0,0.3183098861837,3.0]
end
@btime hh()

function no_prealloc(A,B)
    C = rand(1000,1000) # preallocate array inside of function
    for j in 1:1000, i in 1:1000
        val = A[i,j] + B[i,j]
        C[i,j] = val
    end
end
@btime no_prealloc(A,B) # 3 ms, 2 allocations

C = zeros(1000,1000) # preallocate array outside of function
function w_prealloc!(C,A,B)
    for j in 1:1000, i in 1:1000
        val = A[i,j] + B[i,j]
        C[i,j] = val
    end
end
@btime w_prealloc!(C,A,B) # 2 ms, 0 allocations

w_prealloc!(C,A,B)

C

a = [0 0 0]
v = [1 2 3]
function f!(a,v)
    for i=1:100
        a .+= v
    end
end

@btime f!(a,v)

v = @SVector [1, 2, 3]
function g!(a,v)
    for i=1:100
        a' .+= v
    end
end

@btime g!(a,v)

A.+B.+C./10

function unfused(A,B,C)
    A + B + C
end
@btime unfused(A,B,C) # same as no_broadcast()
D
function fused(A,B,C,D)
    tmp = A .+ B
     tmp .+= C
end
@btime fused(A,B,C)

b = rand(100)
function dothis(b)
    a = sin.(b)
end

dothis(b)


text:
using LIBTIFF
using Base.Test

# test reading

grey_float32 = fill(Float32(0), (50, 60))
grey_float32[20:40, 10:15] = 1
grey_int = fill(Int8(0), (50, 60))
grey_int[20:40, 10:15] = -1

@test all(tiffread("test_data/grey_tile_float32.tif") .== grey_float32)
@test all(tiffread("test_data/grey_tile_float16.tif") .== Float16.(grey_float32))
@test all(tiffread("test_data/grey_tile_uint16.tif") .== UInt16.(grey_float32) * (2^16 - 1))
@test all(tiffread("test_data/grey_tile_uint32.tif") .== UInt32.(grey_float32) * (2^32 - 1))

@test all(tiffread("test_data/grey_scanline_float16.tif") .== Float16.(grey_float32))
@test all(tiffread("test_data/grey_scanline_uint8.tif") .== UInt8.(grey_float32) * 255)
@test all(tiffread("test_data/grey_scanline_int8.tif") .== grey_int)

@test all(grey_float32[30:48, 7:25] .== tiffopen("test_data/grey_tile_float32.tif") do file
        file[30:48, 7:25]
    end)

@test all(Float16.(grey_float32[17:45, 2:56]) .== tiffopen("test_data/grey_scanline_float16.tif") do file
              file[17:45, 2:56]
          end)


text:
module InvariantMeasures

using IntervalArithmetic: range_atan
abstract type NormKind end
struct L1 <: NormKind end
struct L2 <: NormKind end
struct Linf <: NormKind end
struct Lipschitz <: NormKind end
struct TotalVariation <: NormKind end
struct ℓ1 <: NormKind end
struct ℓinf <: NormKind end

# the module Contractors does not depend on any submodule
include("Contractors.jl")

include("DynamicDefinition.jl")
include("BasisDefinition.jl")
include("GenericAssembler.jl")
include("GenericEstimate.jl")

include("PwDynamicDefinition.jl")
include("Mod1Dynamic.jl")
include("Mod1PwDynamic.jl")
include("IterateDynamic.jl")

include("UlamBasis.jl")
include("HatBasis.jl")

using .DynamicDefinition, .BasisDefinition, .Mod1DynamicDefinition, .Contractors, .PwDynamicDefinition

include("Norms.jl")
include("pitrig.jl")
include("NormsOfPowers.jl")

include("preimages.jl")

include("precompile.jl")

export NormKind, L1, Linf, Lipschitz, TotalVariation

export Ulam, PwMap, Mod1Dynamic, Basis, Dynamic, assemble, preim, Hat,
	EquispacedPartition, norms_of_powers, sinpi, cospi, dfly,
	DiscretizedOperator, IntegralPreservingDiscretizedOperator, NonIntegralPreservingDiscretizedOperator,
	opnormbound, weak_norm, strong_norm, aux_norm, integral_covector, distance_from_invariant,
	mod1_dynamic, Iterate, derivative, distorsion, endpoints, nbranches, branch,
	expansivity, max_distorsion, is_refinement,
	skip_beginning, last_end, preimages

import ValidatedNumerics: Interval
export Interval

include("C2Basis.jl")
using .C2BasisDefinition
export C2Basis
include("ContractionC1.jl")

include("ConvergenceRatesOriginal.jl")
export convergencerateabstract

include("HigherDFLY.jl")
include("Chebyshev.jl")

# a special example, the induced map for the LSV map
include("InducedLSV.jl")
using .InducedLSVMapDefinition
export ApproxInducedLSV

include("NoiseKernel.jl")
export UniformNoiseUlam
include("NormsOfPowersNoise.jl")
export powernormboundsnoise, finepowernormboundsnoise, abstractpowernormboundsnoise, invariant_vector_noise, distance_from_invariant_noise

include("Observables.jl")
export Observable, discretizationlogder, integrateobservable

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# 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 holder 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.
# *********************************************************************************
"""
    add_prod_incent_vars_and_constraints(m, p)

When pbi_techs is not empty this function is called to add the variables and constraints for modeling production based
incentives.
"""
function add_prod_incent_vars_and_constraints(m, p)

    @variable(m, dvProdIncent[p.techs.pbi] >= 0)
    @variable(m, binProdIncent[p.techs.pbi], Bin)

	##Constraint (6a)-1: Production Incentive Upper Bound
	@constraint(m, ProdIncentUBCon[t in p.techs.pbi],
		m[:dvProdIncent][t] <= m[:binProdIncent][t] * p.pbi_max_benefit[t] * p.pbi_pwf[t] * p.third_party_factor)

	##Constraint (6a)-2: Production Incentive According to Production
	@constraint(m, IncentByProductionCon[t in p.techs.pbi],
		m[:dvProdIncent][t] <= p.hours_per_time_step * p.pbi_benefit_per_kwh[t] * p.pbi_pwf[t] * p.third_party_factor *
			sum(p.production_factor[t, ts] * m[:dvRatedProduction][t,ts] for ts in p.time_steps)
	)
	##Constraint (6b): System size max to achieve production incentive
	@constraint(m, IncentBySystemSizeCon[t in p.techs.pbi],
		m[:dvSize][t]  <= p.pbi_max_kw[t] + p.max_sizes[t] * (1 - m[:binProdIncent][t])
    )

	m[:TotalProductionIncentive] = @expression(m, sum(m[:dvProdIncent][t] for t in p.techs.pbi))
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
module GradientSketching

using LinearAlgebra, Random, CatViews

export BiasSEGA, SEGA
export project!, projecta!, gradient, unbias!

### General projection methods

"""
    project!(h::AbstractVector, s∇::Number, s::AbstractVector; Binv=I)

Project `h` onto the space consisting of all `y` such that `s'*y = s∇`. Updates `h` in-place.

"""
function project!(h::AbstractVector, s∇::Number, s::AbstractVector; Binv=I)
    length(s) == length(h) || throw(DimensionMismatch("s has length $(length(s)), h has length $(length(h))"))
    Binv == I || size(Binv) == (length(s), length(s)) || throw(DimensionMismatch("Binv has dimensions $(size(Binv)), s has dimensions ($(length(s)),1)"))
    P = Binv * (s ./ dot(s, Binv, s))
    c = dot(s, h) - s∇
    for i in 1:length(h)
        h[i] -= c*P[i]
    end
    h
end

function project!(h::AbstractMatrix, s∇::AbstractVector, s::AbstractVector; Binv=I) where T
    length(s) == size(h, 1) || throw(DimensionMismatch("s has length $(length(s)), h has dimensions $(size(h))"))
    length(s∇) == size(h, 2) || throw(DimensionMismatch("s∇ has length $(length(s∇)), h has dimensions $(size(h))"))
    Binv == I || size(Binv) == (length(s), length(s)) || throw(DimensionMismatch("Binv has dimensions $(size(Binv)), s has dimensions ($(length(s)),1)"))
    P = Binv * (s ./ dot(s, Binv, s))
    for col in 1:size(h, 2)
        hv = view(h, :, col)
        c  = dot(s, hv) - s∇[col]
        for row in 1:length(hv)
            hv[row] -= c*P[row]
        end
    end
    h
end

function project!(h::AbstractVector, S∇::AbstractVector, S::AbstractMatrix; Binv=I)
    size(S, 1) == size(h, 1) || throw(DimensionMismatch("S has dimensions $(size(S)), h has dimensions $(size(h))"))
    size(S∇, 2) == size(h, 2) || throw(DimensionMismatch("S∇ has dimensions $(size(S∇)), h has dimensions $(size(h))"))
    Binv == I || error("not implemented")
    Binv == I || size(Binv) == (size(S, 1), size(S, 1)) || throw(DimensionMismatch("Binv has dimensions $(size(Binv)), S has dimensions $(size(S))"))
    h .-= S * (S \ (S' \ (S'*h .- S∇)))
end

function project!(h::AbstractMatrix, S∇::AbstractMatrix, S::AbstractMatrix; Binv=I)
    size(S, 1) == size(h, 1) || throw(DimensionMismatch("S has dimensions $(size(S)), h has dimensions $(size(h))"))
    size(S∇, 2) == size(h, 2) || throw(DimensionMismatch("S∇ has dimensions $(size(S∇)), h has dimensions $(size(h))"))
    Binv == I || error("not implemented")
    Binv == I || size(Binv) == (size(S, 1), size(S, 1)) || throw(DimensionMismatch("Binv has dimensions $(size(Binv)), S has dimensions $(size(S))"))
    h .-= S * (S \ (S' \ (S'*h .- S∇)))
end

"""
    project!(h::AbstractVector, i::Integer, S)

Coordinate sketch projection.
"""
project!(h::AbstractVector, S∇::Number, i::Integer) = h[i] = S∇

project!(h::AbstractVector{<:AbstractArray}, S∇::AbstractArray, i::Integer) = h[i] .= S∇

function project!(h::AbstractVector{T1}, S∇::AbstractVector{T2}, S::AbstractMatrix; Binv=I) where {T1<:AbstractArray{Tv1,N},T2<:AbstractArray{Tv2,N}} where {Tv1,Tv2,N}
    project!(
        reshape(CatView(h...), length(h[1]), length(h))', # convert h to a matrix by unrolling the component arrays
        reshape(CatView(S∇...), length(S∇[1]), length(S∇))',
        S, Binv=Binv,
    )
    h
end

function project!(h::AbstractVector{T1}, S∇::T2, S::AbstractVector; Binv=I) where {T1<:AbstractArray{Tv1,N},T2<:AbstractArray{Tv2,N}} where {Tv1,Tv2,N}
    project!(
        reshape(CatView(h...), length(h[1]), length(h))', # convert h to a matrix by unrolling the component arrays
        reshape(S∇, length(S∇)), # convert to a vector
        S, Binv=Binv,
    )
    h
end

"""
    projecta!(h::AbstractArray, S∇, S::AbstractMatrix; Binv=I, γ::Integer=5)

Approximate version of `project!` that projects onto each column of `S` separately. Iterates over a random permutation
of the columns `γ` times. This method is equivalent to `project!` if the columns of `S` are orthogonal.
"""
function projecta!(h::AbstractArray, S∇, S::AbstractMatrix; Binv=I, γ::Integer=5)
    γ >= 1 || throw(DomainError(γ, "γ must be positive."))
    p = collect(1:size(S, 1))
    for _ in 1:γ
        shuffle!(p)
        for i in p
            project!(h, selectdim(S∇, 1, i), view(S, :, i), Binv=Binv)
        end
    end
    h
end

"""
    projecta!(h::AbstractArray, S∇::AbstractVector, S::AbstractMatrix; Binv=I, γ::Integer=5)

Approximate version of `project!` for vector gradients.
"""
function projecta!(h::AbstractArray, S∇::AbstractVector, S::AbstractMatrix; Binv=I, γ::Integer=5)
    γ >= 1 || throw(DomainError(γ, "must be positive."))
    p = collect(1:size(S, 1))
    for _ in 1:γ
        shuffle!(p)
        for i in p
            project!(h, S∇[i], view(S, :, i), Binv=Binv)
        end
    end
    h
end

projecta!(h::AbstractArray, S∇, s::AbstractVector; Binv=I, γ::Integer=5) = project!(h, S∇, s; Binv=Binv)

### Biased SEGA

"""

Biased SEGA gradient estimator

"""
struct BiasSEGA{T,N}
    h::Array{T,N}
end

BiasSEGA{T}(dims::Tuple) where T = BiasSEGA{T,length(dims)}(zeros(T, dims))
BiasSEGA(dims::Tuple) = BiasSEGA{Float64}(dims)
BiasSEGA{T}(dim::Integer) where T = BiasSEGA{T}((dim,))
BiasSEGA(dim::Integer) = BiasSEGA{Float64}((dim,))

Base.eltype(sega::BiasSEGA{T}) where T = T
Base.size(sega::BiasSEGA) = size(sega.h)
Base.size(sega::BiasSEGA, i::Integer) = size(sega.h, i)
Base.show(io::IO, sega::BiasSEGA) = print(io, "BiasSEGA{$(eltype(sega)),$(size(sega))")

"""
    project!(sega::BiasSEGA, args...; kwargs...)

Calls `project!(sega.h, args...; kwargs...)`.

"""
project!(sega::BiasSEGA, args...; kwargs...) = project!(sega.h, args...; kwargs...)
projecta!(sega::BiasSEGA, args...; kwargs...) = projecta!(sega.h, args...; kwargs...)

gradient(sega::BiasSEGA) = sega.h

### Unbiased SEGA

"""

SEGA gradient estimator

"""
mutable struct SEGA{T,N}
    τ::Float64          # = 1/θ, where θ is the un-biasing coefficient
    hp::Array{T,N}      # Previous estimate
    h::Array{T,N}       # Biased gradient estimate
    g::Array{T,N}       # Unbiased gradient estimate
end

SEGA{T}(dims::Tuple) where T = SEGA(0.0, zeros(T, dims), zeros(T, dims), zeros(T, dims))
SEGA(dims::Tuple) = SEGA{Float64}(dims)
SEGA{T}(dim::Integer) where T = SEGA{T}((dim,))
SEGA(dim::Integer) = SEGA{Float64}((dim,))

function SEGA(h::Array{<:Number}) where {T<:Number}
    hp = copy(h)
    g = similar(h)
    g .= 0
    SEGA(0.0, hp, h, g)
end

function SEGA(h::Array{<:AbstractArray,N}) where {T,N}
    hp = [copy(hi) for hi in h]
    g = [similar(hi) for hi in h]
    for gi in g
        gi .= 0
    end
    SEGA(0.0, hp, h, g)
end

function SEGA(hp::Array{T,N}, h::Array{T,N}, g::Array{T,N}) where {T,N}
    size(hp) == size(h) || throw(DimensionMismatch("hp has dimensions $(size(hp)), h has dimensions $(size(h))"))
    size(g) == size(h) || throw(DimensionMismatch("g has dimensions $(size(g)), h has dimensions $(size(h))"))
    SEGA(0.0, hp, h, g)
end

Base.eltype(sega::SEGA{T}) where T = T
Base.size(sega::SEGA) = size(sega.h)
Base.size(sega::SEGA, i::Integer) = size(sega.h, i)
Base.show(io::IO, sega::SEGA) = print(io, "SEGA{$(eltype(sega)),$(size(sega))")

function project!(sega::SEGA, τ::Real, args...; kwargs...)
    sega.τ += τ
    project!(sega.h, args...; kwargs...)
end

function projecta!(sega::SEGA, θ::Real, args...; kwargs...)
    sega.τ += τ
    projecta!(sega.h, args...; kwargs...)
end

"""
    unbias!(sega::SEGA)

Compute an unbiased estimate of the gradient.
"""
function unbias!(sega::SEGA)
    0 < sega.τ || throw(DomainError(sega.τ, "τ must be positive"))
    θ = 1/sega.τ
    sega.g .= (1-θ).*sega.hp .+ θ.*sega.h
    sega.hp .= sega.h
    sega.τ = 0.0
    sega.g
end

function unbias!(sega::SEGA{<:AbstractArray})
    0 < sega.τ || throw(DomainError(sega.τ, "τ must be positive"))
    θ = 1/sega.τ
    for (gi, hpi, hi) in zip(sega.g, sega.hp, sega.h)
        gi .= (1-θ).*hpi .+ θ.*hi
        hpi .= hi
    end
    sega.τ = 0.0
    sega.g
end

gradient(sega::SEGA) = sega.g

end


text:
# should always be specified for custom array/tensor types
function similarstructure_from_indices end

# generic definition, net very efficient, provide more efficient version if possible
memsize(a::Any) = Base.summarysize(a)

# generic definitions, should be overwritten if your array/tensor type does not support
# Base.similar(object, eltype, structure)
function similar_from_indices(T::Type, p1::IndexTuple, p2::IndexTuple, A, CA::Symbol)
    structure = similarstructure_from_indices(T, p1, p2, A, CA)
    similar(A, T, structure)
end
function similar_from_indices(T::Type, poA::IndexTuple, poB::IndexTuple,
                                p1::IndexTuple, p2::IndexTuple,
                                A, B, CA::Symbol, CB::Symbol)
    structure = similarstructure_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
    similar(A, T, structure)
end

# should work generically but can be overwritten
Base.@pure function similartype_from_indices(T::Type, p1, p2, A, CA)
    Core.Compiler.return_type(similar_from_indices,
                                Tuple{Type{T}, typeof(p1), typeof(p2), typeof(A), Symbol})
end
Base.@pure function similartype_from_indices(T::Type, poA, poB, p1, p2, A, B, CA, CB)
    Core.Compiler.return_type(similar_from_indices,
                                Tuple{Type{T}, typeof(poA), typeof(poB),
                                        typeof(p1), typeof(p2), typeof(A), typeof(B),
                                        Symbol, Symbol})
end

# generic, should probably not be overwritten
function cached_similar_from_indices(sym::Symbol, T::Type,
                                        p1::IndexTuple, p2::IndexTuple,
                                        A, CA::Symbol)
    if use_cache()
        structure = similarstructure_from_indices(T, p1, p2, A, CA)
        typ = similartype_from_indices(T, p1, p2, A, CA)
        key = (sym, taskid(), typ, structure)
        C::typ = get!(cache, key) do
            similar_from_indices(T, p1, p2, A, CA)
        end
        return C
    else
        return similar_from_indices(T, p1, p2, A, CA)
    end
end
function cached_similar_from_indices(sym::Symbol, T::Type,
                                        poA::IndexTuple, poB::IndexTuple,
                                        p1::IndexTuple, p2::IndexTuple,
                                        A, B, CA::Symbol, CB::Symbol)

    if use_cache()
        structure = similarstructure_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
        typ = similartype_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
        key = (sym, taskid(), typ, structure)
        C::typ = get!(cache, key) do
            similar_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
        end
        return C
    else
        return similar_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
    end
end


text:
relu(x::AbstractArray) = max.(x, 0)
relu_grad(x::AbstractArray) = float(x .> 0)

@primitive relu(x::AbstractArray)
@primitive relu_grad(x::AbstractArray)

@grad relu(x::AbstractArray) 1 relu_grad(x)

relu_test(x) = sum(relu(x))

@testset "macro" begin

    tape = Tape()
    x = tracked(tape, rand(10))
    y = relu(x)

    @test length(tape) == 1
    @test tape[1].val == x.val

    grad(relu_test, rand(5))

end


text:
function birkdecomp(X, ε=1e-12; max_rep = 1)

    n = size(X,1);                                  # get size of Birkhoff polytope
    x_star = reshape(X, n*n);                       # reshape doubly stochastic to vector
    B = birkhoffPolytope(n);                        # Birkhoff polytope
    ε = max(ε,1e-15);                         # fix the maximum minimum precision
    max_iter = (n-1)^2 + 1;

    x = zeros(n*n);                                 # initial point

    extreme_points = zeros(n*n,max_iter);           # extreme points matrix
    θ = zeros(max_iter);                           # weights vector
    approx = Inf;                                   # approximation error
    i = 1;                                          # iteration index

    # Algorithms
    while(approx > ε)

        y = getEPBplus(x_star, x, B, max_rep, ε)
        θi = getBirkhoffStepSize(x_star,x,y)
        x = x + θi*y;
        θ[i] = θi;

        approx = sqrt(sum((abs.(x_star-x)).^2));
        extreme_points[:,i] = y;
        i = i + 1;

    end

    return extreme_points[:,1:i-1], θ[1:i-1]

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SpatialPredicatePartitioner(pred)

A method for partitioning spatial objects with a given spatial
predicate function `pred(x, y)`.
"""
struct SpatialPredicatePartitioner <: AbstractSpatialPredicatePartitioner
  pred::Function
end

(p::SpatialPredicatePartitioner)(x, y) = p.pred(x, y)


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

include("strings/substring.jl")
include("strings/basic.jl")
include("strings/search.jl")
include("strings/unicode.jl")
include("strings/util.jl")
include("strings/io.jl")

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
"""
abstract type TM1Type end

A `TM1Type` is a Julia type representation of a JSON object defined by the TM1
API. Generally:

- The fields of these types should correspond to keys in the JSON object. In the event
the JSON object has a "type" key, the corresponding field name used should be `typ`
(since `type` is a reserved word in Julia).

- The method `namefield` should be defined on every `TM1Type`. This method returns the
type's identity in the form used for URI construction.

- A TM1Type's field types should be Union{Nothing, T} of either: concrete types, a
Vectors of concrete types, or Dicts.

"""
abstract type TM1Type end

"""
@tm1def typeexpr

Define a new `TM1Type` specified by `typeexpr`, adding default constructors for
conversions from `Dict`s and keyword arguments.
"""
macro tm1def(expr)
  # a very simplified form of Base.@kwdef
  expr = macroexpand(__module__, expr) # to expand @static
  expr isa Expr && expr.head == :struct && expr.args[2] isa Symbol ||
    error("Invalid usage of @tm1type")
  T = expr.args[2]
  expr.args[2] = :($T <: TM1Type)

  params_ex = Expr(:parameters)
  call_args = Any[]

  for ei in expr.args[3].args
    if ei isa Expr && ei.head == :(::)
      var = ei.args[1]
      S = ei.args[2]
      push!(params_ex.args, Expr(:kw, var, nothing))
      push!(
        call_args,
        :($var === nothing ? $var : prune_tm1_value($var, unwrap_union_types($S))),
      )
    end
  end
  quote
    Base.@__doc__($(esc(expr)))
    ($(esc(T)))($params_ex) = ($(esc(T)))($(call_args...))
    $(esc(T))(data::Dict) = json2tm1($T, data)
  end
end

function Base.:(==)(a::TM1Type, b::TM1Type)
  if typeof(a) != typeof(b)
    return false
  end

  for field in fieldnames(typeof(a))
    aval, bval = getfield(a, field), getfield(b, field)
    if (aval === nothing) == (bval === nothing)
      if aval !== nothing && aval != bval
        return false
      end
    else
      return false
    end
  end

  return true
end

# `namefield` is overloaded by various TM1Types to allow for more generic
# input to AP functions that require a name to construct URI paths via `name`
name(val) = val
name(t::TM1Type) = namefield(t)

########################################
# Converting JSON Dicts to TM1Types #
########################################

# Unwrap Union{Nothing, Foo} to just Foo
unwrap_union_types(T) = T
function unwrap_union_types(T::Union)
  if T.a == Nothing
    return T.b
  end
  return T.a
end

function extract_nullable(data::Dict, key, ::Type{T}) where {T}
  if haskey(data, key)
    val = data[key]
    if val !== nothing
      if T <: Vector
        V = eltype(T)
        return V[prune_tm1_value(v, unwrap_union_types(V)) for v in val]
      else
        return prune_tm1_value(val, unwrap_union_types(T))
      end
    end
  end
  return nothing
end

prune_tm1_value(val::T, ::Type{Any}) where {T} = val
prune_tm1_value(val::T, ::Type{T}) where {T} = val
prune_tm1_value(val, ::Type{T}) where {T} = T(val)
prune_tm1_value(val::AbstractString, ::Type{Dates.DateTime}) = Dates.DateTime(chopz(val))

# ISO 8601 allows for a trailing 'Z' to indicate that the given time is UTC.
# Julia's Dates.DateTime constructor doesn't support this, but TM1's time
# strings may (??) contain it. This method ensures that a string's trailing 'Z',
# if present, has been removed.

# Not sure this is necessary
function chopz(str::AbstractString)
  if !(isempty(str)) && last(str) == 'Z'
    return chop(str)
  end
  return str
end

# Calling `json2tm1(::Type{TM<:TM1Type}, data::Dict)` will parse the given
# dictionary into the type `TM` with the expectation that the fieldnames of
# `TM` are keys of `data`, and the corresponding values can be converted to the
# given field types.
@generated function json2tm1(::Type{TM}, data::Dict) where {TM<:TM1Type}
  types = unwrap_union_types.(collect(TM.types))
  fields = fieldnames(TM)
  args = Vector{Expr}(undef, length(fields))
  for i in eachindex(fields)
    field, T = fields[i], types[i]
    key = field == :typ ? "type" : string(field)
    args[i] = :(extract_nullable(data, $key, $T))
  end
  return :(TM($(args...))::TM)
end

#############################################
# Converting TM1Type Dicts to JSON Dicts #
#############################################

tm12json(val) = val
tm12json(uri::HTTP.URI) = string(uri)
tm12json(dt::Dates.DateTime) = string(dt) * "Z"
tm12json(v::Vector) = [tm12json(i) for i in v]

function tm12json(t::TM1Type)
  results = Dict()
  for field in fieldnames(typeof(t))
    val = getfield(t, field)
    if val !== nothing
      key = field == :typ ? "type" : string(field)
      results[key] = tm12json(val)
    end
  end
  return results
end

function tm12json(data::Dict{K}) where {K}
  results = Dict{K,Any}()
  for (key, val) in data
    results[key] = tm12json(val)
  end
  return results
end

###################
# Pretty Printing #
###################

function Base.show(io::IO, t::TM1Type)
  if get(io, :compact, false)
    uri_id = namefield(t)
    if uri_id === nothing
      print(io, typeof(t), "(…)")
    else
      print(io, typeof(t), "($(repr(uri_id)))")
    end
  else
    print(io, "$(typeof(t)) (all fields are Union{Nothing, T}):")
    for field in fieldnames(typeof(t))
      val = getfield(t, field)
      if !(val === nothing)
        println(io)
        print(io, "  $field: ")
        if isa(val, Vector)
          print(io, typeof(val))
        else
          show(IOContext(io, :compact => true), val)
        end
      end
    end
  end
end


text:
@testset "Fields" begin
    @test isa(ℝ, Field)
    @test isa(ℂ, Field)

    for T in (Int8, Int16, Int32, Int64, BigInt)
        @test one(T) ∈ ℝ
        @test one(Rational{T}) ∈ ℝ
        @test !(one(Complex{T}) ∈ ℝ)
        @test !(one(Complex{Rational{T}}) ∈ ℝ)
        @test one(T) ∈ ℂ
        @test one(Rational{T}) ∈ ℂ
        @test one(Complex{T}) ∈ ℂ
        @test one(Complex{Rational{T}} ∈ ℂ)

        @test T ⊆ ℝ
        @test Rational{T} ⊆ ℝ
        @test !(Complex{T} ⊆ ℝ)
        @test !(Complex{Rational{T}} ⊆ ℝ)
        @test T ⊆ ℂ
        @test Rational{T} ⊆ ℂ
        @test Complex{T} ⊆ ℂ
        @test Complex{Rational{T}} ⊆ ℂ
    end
    for T in (Float32, Float64, BigFloat)
        @test one(T) ∈ ℝ
        @test !(one(Complex{T}) ∈ ℝ)
        @test one(T) ∈ ℂ
        @test one(Complex{T} ∈ ℂ)

        @test T ⊆ ℝ
        @test !(Complex{T} ⊆ ℝ)
        @test T ⊆ ℂ
        @test Complex{T} ⊆ ℂ
    end
end

@testset "ElementarySpace: CartesianSpace" begin
    d = 2
    V = ℝ^d
    @test eval(Meta.parse(sprint(show,V))) == V
    @test eval(Meta.parse(sprint(show,typeof(V)))) == typeof(V)
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, CartesianSpace)
    @test V == @inferred(dual(V)) == @inferred(conj(V)) == @inferred(adjoint(V))
    @test field(V) == ℝ
    @test @inferred(sectortype(V)) == Trivial
    @test @inferred(TensorKit.hassector(V, Trivial()))
    @test @inferred(dim(V)) == d == @inferred(dim(V, Trivial()))
    @test @inferred(TensorKit.axes(V)) == Base.OneTo(d)
    @test V == ℝ[d] == ℝ[](d) == typeof(V)(d)
    @test ⊕(V,V) == ℝ^(2d)
    @test @inferred min(V, ℝ^3) == V
    @test @inferred max(V', ℝ^3) == ℝ^3
end

@testset "ElementarySpace: ComplexSpace" begin
    d = 2
    V = ℂ^d
    @test eval(Meta.parse(sprint(show,V))) == V
    @test eval(Meta.parse(sprint(show,typeof(V)))) == typeof(V)
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, ComplexSpace)
    @test @inferred(dual(V)) == @inferred(conj(V)) == @inferred(adjoint(V)) != V
    @test @inferred(field(V)) == ℂ
    @test @inferred(sectortype(V)) == Trivial
    @test @inferred(TensorKit.hassector(V, Trivial()))
    @test @inferred(dim(V)) == d == @inferred(dim(V, Trivial()))
    @test @inferred(TensorKit.axes(V)) == Base.OneTo(d)
    @test V == ℂ[d] == ℂ[](d) == typeof(V)(d)
    @test @inferred(⊕(V,V)) == ℂ^(2d)
    @test_throws SpaceMismatch (⊕(V, V'))
    @test_throws MethodError (⊕(ℝ^d, ℂ^d))
    @test @inferred min(V, ℂ^3) == V
    @test @inferred max(V', ℂ[3]') == ℂ[3]'
end

@testset "ElementarySpace: GeneralSpace" begin
    d = 2
    V = GeneralSpace{ℂ}(d)
    @test eval(Meta.parse(sprint(show,V))) == V
    @test eval(Meta.parse(sprint(show,typeof(V)))) == typeof(V)
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test !isa(V, InnerProductSpace)
    @test !isa(V, EuclideanSpace)
    @test @inferred(dual(V)) != @inferred(conj(V)) != V
    @test @inferred(field(V)) == ℂ
    @test @inferred(sectortype(V)) == Trivial
    @test @inferred(TensorKit.hassector(V, Trivial()))
    @test @inferred(dim(V)) == d == @inferred(dim(V, Trivial()))
    @test @inferred(TensorKit.axes(V)) == Base.OneTo(d)
end

@testset "ElementarySpace: RepresentationSpace{$G}" for G in (ℤ₂, ℤ₃, ℤ₄, U₁, CU₁, SU₂, FibonacciAnyon, ℤ₃ × ℤ₄, U₁ × SU₂, SU₂ × SU₂, ℤ₂ × FibonacciAnyon × FibonacciAnyon)
    if Base.IteratorSize(values(G)) === Base.IsInfinite()
        set = unique(vcat(one(G), [randsector(G) for k = 1:10]))
        gen = (c=>2 for c in set)
    else
        gen = (values(G)[k]=>(k+1) for k in 1:length(values(G)))
    end
    V = RepresentationSpace(gen)
    @test eval(Meta.parse(sprint(show,V))) == V
    @test V' == RepresentationSpace(gen; dual = true)
    @test V == @inferred RepresentationSpace(gen...)
    @test V' == @inferred RepresentationSpace(gen...; dual = true)
    @test V == @inferred RepresentationSpace(tuple(gen...))
    @test V' == @inferred RepresentationSpace(tuple(gen...); dual = true)
    @test V == @inferred RepresentationSpace{G}(gen)
    @test V' == @inferred RepresentationSpace{G}(gen; dual = true)
    @test V == @inferred RepresentationSpace{G}(gen...)
    @test V' == @inferred RepresentationSpace{G}(gen...; dual = true)
    @test eval(Meta.parse(sprint(show,V))) == V
    @test eval(Meta.parse(sprint(show,typeof(V)))) == typeof(V)
    W = RepresentationSpace(one(G)=>1) # space with a single sector
    @test W == RepresentationSpace(one(G)=>1, randsector(G) => 0)
    # randsector never returns trivial sector, so this cannot error
    @test_throws ArgumentError RepresentationSpace(one(G)=>1, randsector(G) => 0, one(G)=>3)
    @test eval(Meta.parse(sprint(show,W))) == W
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, RepresentationSpace)
    @test isa(V, RepresentationSpace{G})
    @test isa(V, Base.IteratorSize(values(G)) == Base.IsInfinite() ?
                    TensorKit.GenericRepresentationSpace{G} :
                    TensorKit.FiniteRepresentationSpace{G})
    @test @inferred(dual(V)) == @inferred(conj(V)) == @inferred(adjoint(V)) != V
    @test @inferred(field(V)) == ℂ
    @test @inferred(sectortype(V)) == G
    slist = @inferred sectors(V)
    @test @inferred(TensorKit.hassector(V, first(slist)))
    @test @inferred(dim(V)) == sum((@inferred(dim(s)*dim(V,s))) for s in slist)
    if hasfusiontensor(G)
        @test @inferred(TensorKit.axes(V)) == Base.OneTo(dim(V))
    end
    @test @inferred(⊕(V,V)) == @inferred max(V, ⊕(V,V))
    @test V == @inferred min(V, ⊕(V,V))
    @test min(V, RepresentationSpace(one(G)=>3)) == RepresentationSpace(one(G)=>2)
    @test_throws SpaceMismatch (⊕(V, V'))
end