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

##############################################
#
#   Time Domain Reduction
#    - Jack Morris
#
# Use kmeans or kemoids to cluster raw load profiles and resource capacity factor profiles
# into representative periods. Use Extreme Periods to capture noteworthy periods or
# periods with notably poor fits.
#
#  Inputs
#
#  In time_domain_reduction_settings.yml, include the following:
#
#  - Timesteps_per_Rep_Period - Typically 168 timesteps (e.g., hours) per period, this designates the length
#     of each representative period.
#  - UseExtremePeriods - Either 1 or 0, this designates whether or not to include
#     outliers (by performance or load/resource extreme) as their own representative periods.
#     This setting automatically includes the periods with maximum load, minimum solar cf and
#     minimum wind cf as extreme periods.
#  - ClusterMethod - Either 'kmeans' or 'kmedoids', this designates the method used to cluster
#     periods and determine each point's representative period.
#  - ScalingMethod - Either 'N' or 'S', this designates directs the module to normalize ([0,1])
#     or standardize (mean 0, variance 1) the input data.
#  - MinPeriods - The minimum number of periods used to represent the input data. If using
#     UseExtremePeriods, this must be at least three. If IterativelyAddPeriods if off,
#    this will be the total number of periods.
#  - MaxPeriods - The maximum number of periods - both clustered periods and extreme periods -
#     that may be used to represent the input data.
#  - IterativelyAddPeriods - Either 1 or 0, this designates whether or not to add periods
#     until the error threshold between input data and represented data is met or the maximum
#     number of periods is reached.
#  - Threshold - Iterative period addition will end if the period farthest (Euclidean Distance)
#     from its representative period is within this percentage of the total possible error (for normalization)
#     or ~95% of the total possible error (for standardization). E.g., for a threshold of 0.01,
#     every period must be within 1% of the spread of possible error before the clustering
#     iterations will terminate (or until the max number of periods is reached).
#  - IterateMethod - Either 'cluster' (Default) or 'extreme', this designates whether to add clusters to
#     the kmeans/kmedoids method or to set aside the worst-fitting periods as a new extreme periods.
#  - nReps - Default 200, the number of times to repeat each kmeans/kmedoids clustering at the same setting.
#  - LoadWeight - Default 1, this is an optional multiplier on load columns in order to prioritize
#     better fits for load profiles over resource capacity factor profiles.
#  - WeightTotal - Default 8760, the sum to which the relative weights of representative periods will be scaled.
#  - ClusterFuelPrices - Either 1 or 0, this indicates whether or not to use the fuel price
#     time series in Fuels_data.csv in the clustering process. If 'no', this function will still write
#     Fuels_data_clustered.csv with reshaped fuel prices based on the number and size of the
#     representative weeks, assuming a constant time series of fuel prices with length equal to the
#     number of timesteps in the raw input data.
#
#
#############################################

# Store the paths of the current working directory, GenX, and Settings YAML file
if isfile("PreCluster.jl")
    genx_path = cd(pwd, "../..") # pwd() grandparent <-- TDR called from PreCluster.jl
else
    genx_path = pwd()         # <-- TDR called from Run_test.jl
end
settings_path = joinpath(genx_path, "GenX_settings.yml")

# Load GenX modules
push!(LOAD_PATH, genx_path)
push!(LOAD_PATH, pwd())

using YAML
using DataFrames
using StatsBase
using Clustering
using Distances
using Documenter
using CSV

@doc raw"""
    rmse_score(y_true, y_pred)

Calculates Root Mean Square Error.

```math
RMSE = \sqrt{\frac{1}{n}\Sigma_{i=1}^{n}{\Big(\frac{d_i -f_i}{\sigma_i}\Big)^2}}
```

"""
function rmse_score(y_true, y_pred)
    errors = y_pred - y_true
    errors² = errors .^ 2
    mse = mean(errors²)
    rmse = sqrt(mse)
    return rmse
end

@doc raw"""
    parse_data(myinputs)

Get load, solar, wind, and other curves from the input data.

"""
function parse_data(myinputs)
    RESOURCES = myinputs["RESOURCE_ZONES"]
    ZONES = myinputs["R_ZONES"]
    # Assuming no missing data
    solar_col_names = []
    wind_col_names = []
    var_col_names = []
    solar_profiles = []
    wind_profiles = []
    var_profiles = []

    # LOAD - Load_data.csv
    load_profiles = [ myinputs["pD"][:,l] for l in 1:size(myinputs["pD"],2) ]
    load_col_names = ["Load_MW_z"*string(l) for l in 1:size(load_profiles)[1]]
    load_zones = [l for l in 1:size(load_profiles)[1]]
    col_to_zone_map = Dict("Load_MW_z"*string(l) => l for l in 1:size(load_profiles)[1])

    # CAPACITY FACTORS - Generators_variability.csv
    for r in 1:length(RESOURCES)
        if occursin("PV", RESOURCES[r]) || occursin("pv", RESOURCES[r]) || occursin("Pv", RESOURCES[r]) || occursin("Solar", RESOURCES[r]) || occursin("SOLAR", RESOURCES[r]) || occursin("solar", RESOURCES[r])
            push!(solar_col_names, RESOURCES[r])
            push!(solar_profiles, myinputs["pP_Max"][r,:])
        elseif occursin("Wind", RESOURCES[r]) || occursin("WIND", RESOURCES[r]) || occursin("wind", RESOURCES[r])
            push!(wind_col_names, RESOURCES[r])
            push!(wind_profiles, myinputs["pP_Max"][r,:])
        end
        push!(var_col_names, RESOURCES[r])
        push!(var_profiles, myinputs["pP_Max"][r,:])
        col_to_zone_map[RESOURCES[r]] = ZONES[r]
    end

    fuel_col_names = string.(myinputs["fuels"])
    fuel_profiles = []
    AllFuelsConst = true
    for f in 1:length(fuel_col_names)
        push!(fuel_profiles, myinputs["fuel_costs"][fuel_col_names[f]])
        if AllFuelsConst && (minimum(myinputs["fuel_costs"][fuel_col_names[f]]) != maximum(myinputs["fuel_costs"][fuel_col_names[f]]))
            AllFuelsConst = false
        end
    end
    all_col_names = [load_col_names; var_col_names; fuel_col_names]
    all_profiles = [load_profiles..., var_profiles..., fuel_profiles...]
    return load_col_names, var_col_names, solar_col_names, wind_col_names, fuel_col_names, all_col_names,
         load_profiles, var_profiles, solar_profiles, wind_profiles, fuel_profiles, all_profiles,
         col_to_zone_map, AllFuelsConst
end

@doc raw"""
    check_condition(Threshold, R, OldColNames, ScalingMethod, TimestepsPerRepPeriod)

Check whether the greatest Euclidean deviation in the input data and the clustered
representation is within a given proportion of the "maximum" possible deviation.

(1 for Normalization covers 100%, 4 for Standardization covers ~95%)

"""
function check_condition(Threshold, R, OldColNames, ScalingMethod, TimestepsPerRepPeriod)
    if ScalingMethod == "N"
        return maximum(R.costs)/(length(OldColNames)*TimestepsPerRepPeriod) < Threshold
    elseif ScalingMethod == "S"
        return maximum(R.costs)/(length(OldColNames)*TimestepsPerRepPeriod*4) < Threshold
    else
        println("INVALID Scaling Method ", ScalingMethod, " / Choose N for Normalization or S for Standardization. Proceeding with N.")
    end
    return maximum(R.costs)/(length(OldColNames)*TimestepsPerRepPeriod) < Threshold
end

@doc raw"""
    get_worst_period_idx(R)

Get the index of the period that is farthest from its representative
period by Euclidean distance.

"""
function get_worst_period_idx(R)
    return argmax(R.costs)
end

@doc raw"""
    cluster(ClusterMethod, ClusteringInputDF, NClusters, nIters)

Get representative periods using cluster centers from kmeans or kmedoids.

K-Means:
https://juliastats.org/Clustering.jl/dev/kmeans.html

K-Medoids:
 https://juliastats.org/Clustering.jl/stable/kmedoids.html
"""
function cluster(ClusterMethod, ClusteringInputDF, NClusters, nIters, v=false)
    if ClusterMethod == "kmeans"
        DistMatrix = pairwise(Euclidean(), Matrix(ClusteringInputDF), dims=2)
        R = kmeans(Matrix(ClusteringInputDF), NClusters, init=:kmcen)

        for i in 1:nIters
            R_i = kmeans(Matrix(ClusteringInputDF), NClusters)

            if R_i.totalcost < R.totalcost
                R = R_i
            end
            if v && (i % (nIters/10) == 0)
                println(string(i) * " : " * string(round(R_i.totalcost, digits=3)) * " " * string(round(R.totalcost, digits=3)) )
            end
        end

        A = R.assignments # get points to clusters mapping - A for Assignments
        W = R.counts # get the cluster sizes - W for Weights
        Centers = R.centers # get the cluster centers - M for Medoids

        M = []
        for i in 1:NClusters
            dists = [euclidean(Centers[:,i], ClusteringInputDF[!, j]) for j in 1:size(ClusteringInputDF, 2)]
            push!(M,argmin(dists))
        end

    elseif ClusterMethod == "kmedoids"
        DistMatrix = pairwise(Euclidean(), Matrix(ClusteringInputDF), dims=2)
        R = kmedoids(DistMatrix, NClusters, init=:kmcen)

        for i in 1:nIters
            R_i = kmedoids(DistMatrix, NClusters)
            if R_i.totalcost < R.totalcost
                R = R_i
            end
            if v && (i % (nIters/10) == 0)
                println(string(i) * " : " * string(round(R_i.totalcost, digits=3)) * " " * string(round(R.totalcost, digits=3)) )
            end
        end

        A = R.assignments # get points to clusters mapping - A for Assignments
        W = R.counts # get the cluster sizes - W for Weights
        M = R.medoids # get the cluster centers - M for Medoids
    else
        println("INVALID ClusterMethod. Select kmeans or kmedoids. Running kmeans instead.")
        return cluster("kmeans", ClusteringInputDF, NClusters, nIters)
    end
    return [R, A, W, M, DistMatrix]
end

@doc raw"""
    RemoveConstCols(all_profiles, all_col_names)

Remove and store the columns that do not vary during the period.

"""
function RemoveConstCols(all_profiles, all_col_names, v=false)
    ConstData = []
    ConstIdx = []
    ConstCols = []
    for c in 1:length(all_col_names)
        Const = minimum(all_profiles[c]) == maximum(all_profiles[c])
        if Const
            if v println("Removing constant col: ", all_col_names[c]) end
            push!(ConstData, all_profiles[c])
            push!(ConstCols, all_col_names[c])
            push!(ConstIdx, c)
        end
    end
    all_profiles = [all_profiles[i] for i in 1:length(all_profiles) if i ∉ ConstIdx]
    all_col_names = [all_col_names[i] for i in 1:length(all_col_names) if i ∉ ConstIdx]
    return all_profiles, all_col_names, ConstData, ConstCols, ConstIdx
end

@doc raw"""

    get_extreme_period(DF, GDF, profKey, typeKey, statKey,
       ConstCols, load_col_names, solar_col_names, wind_col_names)

Identify extreme week by specification of profile type (Load, PV, Wind),
measurement type (absolute (timestep with min/max value) vs. integral
(period with min/max summed value)), and statistic (minimum or maximum).
I.e., the user could want the hour with the most load across the whole
system to be included among the extreme periods. They would select
"Load", "System, "Absolute, and "Max".

"""
function get_extreme_period(DF, GDF, profKey, typeKey, statKey,
    ConstCols, load_col_names, solar_col_names, wind_col_names, v=false)
    if v println(profKey," ", typeKey," ", statKey) end
    if typeKey == "Integral"
        if profKey == "Load"
            (stat, group_idx) = get_integral_extreme(GDF, statKey, load_col_names, ConstCols)
        elseif profKey == "PV"
            (stat, group_idx) = get_integral_extreme(GDF, statKey, solar_col_names, ConstCols)
        elseif profKey == "Wind"
            (stat, group_idx) = get_integral_extreme(GDF, statKey, wind_col_names, ConstCols)
        else
            println("Error: Profile Key ", profKey, " is invalid. Choose `Load', `PV' or `Wind'.")
        end
    elseif typeKey == "Absolute"
        if profKey == "Load"
            (stat, group_idx) = get_absolute_extreme(DF, statKey, load_col_names, ConstCols)
        elseif profKey == "PV"
            (stat, group_idx) = get_absolute_extreme(DF, statKey, solar_col_names, ConstCols)
        elseif profKey == "Wind"
            (stat, group_idx) = get_absolute_extreme(DF, statKey, wind_col_names, ConstCols)
        else
            println("Error: Profile Key ", profKey, " is invalid. Choose `Load', `PV' or `Wind'.")
        end
   else
       println("Error: Type Key ", typeKey, " is invalid. Choose `Absolute' or `Integral'.")
       stat = 0
       group_idx = 0
   end
    return (stat, group_idx)
end

@doc raw"""

    get_integral_extreme(GDF, statKey, col_names, ConstCols)

Get the period index with the minimum or maximum load or capacity factor
summed over the period.

"""
function get_integral_extreme(GDF, statKey, col_names, ConstCols)
    if statKey == "Max"
        (stat, stat_idx) = findmax( sum([GDF[!, Symbol(c)] for c in setdiff(col_names, ConstCols) ]) )
    elseif statKey == "Min"
        (stat, stat_idx) = findmin( sum([GDF[!, Symbol(c)] for c in setdiff(col_names, ConstCols) ]) )
    else
        println("Error: Statistic Key ", statKey, " is invalid. Choose `Max' or `Min'.")
    end
    return (stat, stat_idx)
end

@doc raw"""

    get_absolute_extreme(DF, statKey, col_names, ConstCols)

Get the period index of the single timestep with the minimum or maximum load or capacity factor.

"""
function get_absolute_extreme(DF, statKey, col_names, ConstCols)
    if statKey == "Max"
        (stat, stat_idx) = findmax( sum([DF[!, Symbol(c)] for c in setdiff(col_names, ConstCols) ]) )
        group_idx = DF.Group[stat_idx]
    elseif statKey == "Min"
        (stat, stat_idx) = findmin( sum([DF[!, Symbol(c)] for c in setdiff(col_names, ConstCols) ]) )
        group_idx = DF.Group[stat_idx]
    else
        println("Error: Statistic Key ", statKey, " is invalid. Choose `Max' or `Min'.")
    end
    return (stat, group_idx)
end

@doc raw"""

    scale_weights(W, H)

Linearly scale weights W such that they sum to the desired number of timesteps (hours) H.

```math
w_j \leftarrow H \cdot \frac{w_j}{\sum_i w_i} \: \: \: \forall w_j \in W
```

"""
function scale_weights(W, H, v=false)
    if v println("Weights before scaling: ", W) end
    W = [ float(w)/sum(W) * H for w in W] # Scale to number of hours in input data
    if v
        println("Weights after scaling: ", W)
        println("Sum of Updated Cluster Weights: ", sum(W))
    end
    return W
end

@doc raw"""

    get_load_multipliers(ClusterOutputData, ModifiedData, M, W, LoadCols, TimestepsPerRepPeriod, NewColNames, NClusters, Ncols)

Get multipliers to linearly scale clustered load profiles L zone-wise such that their weighted sum equals the original zonal total load.
Scale load profiles later using these multipliers in order to ensure that a copy of the original load is kept for validation.

Find $k_z$ such that:

```math
\sum_{i \in I} L_{i,z} = \sum_{t \in T, m \in M} C_{t,m,z} \cdot \frac{w_m}{T} \cdot k_z   \: \: \: \forall z \in Z
```

where $Z$ is the set of zones, $I$ is the full time domain, $T$ is the length of one period (e.g., 168 for one week in hours),
$M$ is the set of representative periods, $L_{i,z}$ is the original zonal load profile over time (hour) index $i$, $C_{i,m,z}$ is the
load in timestep $i$ for representative period $m$ in zone $z$, $w_m$ is the weight of the representative period equal to the total number of
hours that one hour in representative period $m$ represents in the original profile, and $k_z$ is the zonal load multiplier returned by the function.

"""
function get_load_multipliers(ClusterOutputData, InputData, M, W, LoadCols, TimestepsPerRepPeriod, NewColNames, NClusters, Ncols, v=false)

    # Compute original zonal total loads
    zone_sums = Dict()
    for loadcol in LoadCols
        zone_sums[loadcol] = sum(InputData[:, loadcol])
    end

    # Compute zonal loads per representative period
    cluster_zone_sums = Dict()
    for m in 1:NClusters
        clustered_lp_DF = DataFrame( Dict( NewColNames[i] => ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i] for i in 1:Ncols if (Symbol(NewColNames[i]) in LoadCols)) )
        cluster_zone_sums[m] = Dict()
        for loadcol in LoadCols
            cluster_zone_sums[m][loadcol] = sum(clustered_lp_DF[:, loadcol])
        end
    end

    # Use representative period weights to compute total zonal load of the representative profile
    # Determine multiplier to bridge the gap between original zonal loads and representative zonal loads
    weighted_cluster_zone_sums = Dict(loadcol => 0.0 for loadcol in LoadCols)
    load_mults = Dict()
    for loadcol in LoadCols
        for m in 1:NClusters
            weighted_cluster_zone_sums[loadcol] += (W[m]/(TimestepsPerRepPeriod))*cluster_zone_sums[m][loadcol]
        end
        load_mults[loadcol] = zone_sums[loadcol]/weighted_cluster_zone_sums[loadcol]
        if v println(loadcol, ": ", weighted_cluster_zone_sums[loadcol], " vs. ", zone_sums[loadcol], " => ", load_mults[loadcol]) end
    end

    # Zone-wise validation that scaled clustered load equals original load (Don't actually scale load in this function)
    if v
        new_zone_sums = Dict(loadcol => 0.0 for loadcol in LoadCols)
        for m in 1:NClusters
            for i in 1:Ncols
                if (NewColNames[i] in LoadCols)
                    # Uncomment this line if we decide to scale load here instead of later. (Also remove "load_mults[NewColNames[i]]*" term from new_zone_sums computation)
                    #ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i] *= load_mults[NewColNames[i]]
                    println("   Scaling ", M[m], " (", NewColNames[i], ") : ", cluster_zone_sums[m][NewColNames[i]], " => ", load_mults[NewColNames[i]]*sum(ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i]))
                    new_zone_sums[NewColNames[i]] += (W[m]/(TimestepsPerRepPeriod))*load_mults[NewColNames[i]]*sum(ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i])
                end
            end
        end
        for loadcol in LoadCols
            println(loadcol, ": ", new_zone_sums[loadcol], " =?= ", zone_sums[loadcol])
        end
    end

    return load_mults
end

@doc raw"""
    cluster_inputs(inpath, settings_path, v=false, norm_plot=false, silh_plot=false, res_plots=false, indiv_plots=false, pair_plots=false)

Use kmeans or kmedoids to cluster raw load profiles and resource capacity factor profiles
into representative periods. Use Extreme Periods to capture noteworthy periods or
periods with notably poor fits.

In Load_data.csv, include the following:

 - Timesteps\_per\_Rep\_Period - Typically 168 timesteps (e.g., hours) per period, this designates the length
     of each representative period.
 - UseExtremePeriods - Either 1 or 0, this designates whether or not to include
    outliers (by performance or load/resource extreme) as their own representative periods.
    This setting automatically includes the periods with maximum load, minimum solar cf and
    minimum wind cf as extreme periods.
 - ClusterMethod - Either 'kmeans' or 'kmedoids', this designates the method used to cluster
    periods and determine each point's representative period.
 - ScalingMethod - Either 'N' or 'S', this designates directs the module to normalize ([0,1])
    or standardize (mean 0, variance 1) the input data.
 - MinPeriods - The minimum number of periods used to represent the input data. If using
    UseExtremePeriods, this must be at least three. If IterativelyAddPeriods if off,
    this will be the total number of periods.
 - MaxPeriods - The maximum number of periods - both clustered periods and extreme periods -
    that may be used to represent the input data.
 - IterativelyAddPeriods - Either 1 or 0, this designates whether or not to add periods
    until the error threshold between input data and represented data is met or the maximum
    number of periods is reached.
 - Threshold - Iterative period addition will end if the period farthest (Euclidean Distance)
    from its representative period is within this percentage of the total possible error (for normalization)
    or ~95% of the total possible error (for standardization). E.g., for a threshold of 0.01,
    every period must be within 1% of the spread of possible error before the clustering
    iterations will terminate (or until the max number of periods is reached).
 - IterateMethod - Either 'cluster' or 'extreme', this designates whether to add clusters to
    the kmeans/kmedoids method or to set aside the worst-fitting periods as a new extreme periods.
 - nReps - The number of times to repeat each kmeans/kmedoids clustering at the same setting.
 - LoadWeight - Default 1, this is an optional multiplier on load columns in order to prioritize
    better fits for load profiles over resource capacity factor profiles.
 - WeightTotal - Default 8760, the sum to which the relative weights of representative periods will be scaled.
 - ClusterFuelPrices - Either 1 or 0, this indicates whether or not to use the fuel price
    time series in Fuels\_data.csv in the clustering process. If 'no', this function will still write
    Fuels\_data\_clustered.csv with reshaped fuel prices based on the number and size of the
    representative weeks, assuming a constant time series of fuel prices with length equal to the
    number of timesteps in the raw input data.
"""
function cluster_inputs(inpath, settings_path, mysetup, v=false)

    if v println(now()) end

    ##### Step 0: Load in settings and data

    # Read time domain reduction settings file time_domain_reduction_settings.yml
    myTDRsetup = YAML.load(open(joinpath(settings_path,"time_domain_reduction_settings.yml")))

    # Accept model parameters from the settings file time_domain_reduction_settings.yml
    TimestepsPerRepPeriod = myTDRsetup["TimestepsPerRepPeriod"]
    ClusterMethod = myTDRsetup["ClusterMethod"]
    ScalingMethod = myTDRsetup["ScalingMethod"]
    MinPeriods = myTDRsetup["MinPeriods"]
    MaxPeriods = myTDRsetup["MaxPeriods"]
    UseExtremePeriods = myTDRsetup["UseExtremePeriods"]
    ExtPeriodSelections = myTDRsetup["ExtremePeriods"]
    Iterate = myTDRsetup["IterativelyAddPeriods"]
    IterateMethod = myTDRsetup["IterateMethod"]
    Threshold = myTDRsetup["Threshold"]
    nReps = myTDRsetup["nReps"]
    LoadWeight = myTDRsetup["LoadWeight"]
    WeightTotal = myTDRsetup["WeightTotal"]
    ClusterFuelPrices = myTDRsetup["ClusterFuelPrices"]
    TimeDomainReductionFolder = mysetup["TimeDomainReductionFolder"]

    # Set output filenames for later
    Load_Outfile = joinpath(TimeDomainReductionFolder, "Load_data.csv")
    GVar_Outfile = joinpath(TimeDomainReductionFolder, "Generators_variability.csv")
    Fuel_Outfile = joinpath(TimeDomainReductionFolder, "Fuels_data.csv")
    PMap_Outfile = joinpath(TimeDomainReductionFolder, "Period_map.csv")
    YAML_Outfile = joinpath(TimeDomainReductionFolder, "time_domain_reduction_settings.yml")

    # Define a local version of the setup so that you can modify the mysetup["ParameterScale"] value to be zero in case it is 1
    mysetup_local = mysetup
    # If ParameterScale =1 then make it zero, since clustered inputs will be scaled prior to generating model
    mysetup_local["ParameterScale"]=0  # Performing cluster and report outputs in user-provided units
    if v println("Loading inputs") end
    myinputs=Dict()
    myinputs = load_inputs(mysetup_local,inpath)
    if v println() end

    # Parse input data into useful structures divided by type (load, wind, solar, fuel, groupings thereof, etc.)
    # TO DO LATER: Replace these with collections of col_names, profiles, zones
    load_col_names, var_col_names, solar_col_names, wind_col_names, fuel_col_names, all_col_names,
         load_profiles, var_profiles, solar_profiles, wind_profiles, fuel_profiles, all_profiles,
         col_to_zone_map, AllFuelsConst = parse_data(myinputs)

    # Remove Constant Columns - Add back later in final output
    all_profiles, all_col_names, ConstData, ConstCols, ConstIdx = RemoveConstCols(all_profiles, all_col_names, v)

    # Determine whether or not to time domain reduce fuel profiles as well based on user choice and file structure (i.e., variable fuels in Fuels_data.csv)
    IncludeFuel = true
    if (ClusterFuelPrices != 1) || (AllFuelsConst) IncludeFuel = false end

    # Put it together!
    InputData = DataFrame( Dict( all_col_names[c]=>all_profiles[c] for c in 1:length(all_col_names) ) )
    if v
        println("Load (MW) and Capacity Factor Profiles: ")
        println(describe(InputData))
        println()
    end
    OldColNames = names(InputData)
    NewColNames = [Symbol.(OldColNames); :GrpWeight]
    Nhours = nrow(InputData) # Timesteps
    Ncols = length(NewColNames) - 1

    ##### Step 1: Normalize or standardize all load, renewables, and fuel data / optionally scale with LoadWeight

    # Normalize/standardize data based on user-provided method
    if ScalingMethod == "N"
        normProfiles = [ StatsBase.transform(fit(UnitRangeTransform, InputData[:,c]; dims=1, unit=true), InputData[:,c]) for c in 1:length(OldColNames)  ]
    elseif ScalingMethod == "S"
        normProfiles = [ StatsBase.transform(fit(ZScoreTransform, InputData[:,c]; dims=1, center=true, scale=true), InputData[:,c]) for c in 1:length(OldColNames)  ]
    else
        println("ERROR InvalidScalingMethod: Use N for Normalization or S for Standardization.")
        println("CONTINUING using 0->1 normalization...")
        normProfiles = [ StatsBase.transform(fit(UnitRangeTransform, InputData[:,c]; dims=1, unit=true), InputData[:,c]) for c in 1:length(OldColNames)  ]
    end

    # Compile newly normalized/standardized profiles
    AnnualTSeriesNormalized = DataFrame(Dict(  OldColNames[c] => normProfiles[c] for c in 1:length(OldColNames) ))

    # Optional pre-scaling of load in order to give it more preference in clutering algorithm
    if LoadWeight != 1   # If we want to value load more/less than capacity factors. Assume nonnegative. LW=1 means no scaling.
        for c in load_col_names
            AnnualTSeriesNormalized[!, Symbol(c)] .= AnnualTSeriesNormalized[!, Symbol(c)] .* LoadWeight
        end
    end

    if v
        println("Load (MW) and Capacity Factor Profiles NORMALIZED! ")
        println(describe(AnnualTSeriesNormalized))
        println()
    end

    ##### STEP 2: Identify extreme periods in the model, Reshape data for clustering

    # Total number of subperiods available in the dataset, where each subperiod length = TimestepsPerRepPeriod
    NumDataPoints = Nhours÷TimestepsPerRepPeriod # 364 weeks in 7 years
    if v println("Total Subperiods in the data set: ", NumDataPoints) end
    InputData[:, :Group] .= (1:Nhours) .÷ (TimestepsPerRepPeriod+0.0001) .+ 1    # Group col identifies the subperiod ID of each hour (e.g., all hours in week 2 have Group=2 if using TimestepsPerRepPeriod=168)

    # Group by period (e.g., week)
    cgdf = combine(groupby(InputData, :Group), [c .=> sum for c in OldColNames])
    cgdf = cgdf[setdiff(1:end, NumDataPoints+1), :]
    rename!(cgdf, [:Group; Symbol.(OldColNames)])

    # Extreme period identification based on user selection in time_domain_reduction_settings.yml
    LoadExtremePeriod = false        # Used when deciding whether or not to scale load curves to equal original total load
    ExtremeWksList = []
    if UseExtremePeriods == 1
      for profKey in keys(ExtPeriodSelections)
          for geoKey in keys(ExtPeriodSelections[profKey])
              for typeKey in keys(ExtPeriodSelections[profKey][geoKey])
                  for statKey in keys(ExtPeriodSelections[profKey][geoKey][typeKey])
                      if ExtPeriodSelections[profKey][geoKey][typeKey][statKey] == 1
                          if profKey == "Load"
                              LoadExtremePeriod = true
                          end
                          if geoKey == "System"
                              if v print(geoKey, " ") end
                              (stat, group_idx) = get_extreme_period(InputData, cgdf, profKey, typeKey, statKey, ConstCols, load_col_names, solar_col_names, wind_col_names, v)
                              push!(ExtremeWksList, floor(Int, group_idx))
                              if v println(group_idx, " : ", stat) end
                          elseif geoKey == "Zone"
                              for z in sort(unique(myinputs["R_ZONES"]))
                                  z_cols = [k for (k,v) in col_to_zone_map if v==z]
                                  if profKey == "Load" z_cols_type = intersect(z_cols, load_col_names)
                                  elseif profKey == "PV" z_cols_type = intersect(z_cols, solar_col_names)
                                  elseif profKey == "Wind" z_cols_type = intersect(z_cols, wind_col_names)
                                  else z_cols_type = []
                                  end
                                  z_cols_type = setdiff(z_cols_type, ConstCols)
                                  if length(z_cols_type) > 0
                                      if v print(geoKey, " ") end
                                      (stat, group_idx) = get_extreme_period(select(InputData, [:Group; Symbol.(z_cols_type)]), select(cgdf, [:Group; Symbol.(z_cols_type)]), profKey, typeKey, statKey, ConstCols, z_cols_type, z_cols_type, z_cols_type, v)
                                      push!(ExtremeWksList, floor(Int, group_idx))
                                      if v println(group_idx, " : ", stat, "(", z, ")") end
                                  else
                                      if v println("Zone ", z, " has no time series profiles of type ", profKey) end
                                  end
                              end
                          else
                              println("Error: Geography Key ", geoKey, " is invalid. Select `System' or `Zone'.")
                          end
                      end
                  end
              end
          end
      end
      if v println(ExtremeWksList) end
      sort!(unique!(ExtremeWksList))
      if v println("Reduced to ", ExtremeWksList) end
    end

    ### DATA MODIFICATION - Shifting InputData and Normalized InputData
    #    from 8760 (# hours) by n (# profiles) DF to
    #    168*n (n period-stacked profiles) by 52 (# periods) DF
    DFsToConcat = [stack(InputData[isequal.(InputData.Group,w),:], OldColNames)[!,:value] for w in 1:NumDataPoints if w <= NumDataPoints ]
    ModifiedData = DataFrame(Dict(Symbol(i) => DFsToConcat[i] for i in 1:NumDataPoints))

    AnnualTSeriesNormalized[:, :Group] .= (1:Nhours) .÷ (TimestepsPerRepPeriod+0.0001) .+ 1
    DFsToConcatNorm = [stack(AnnualTSeriesNormalized[isequal.(AnnualTSeriesNormalized.Group,w),:], OldColNames)[!,:value] for w in 1:NumDataPoints if w <= NumDataPoints ]
    ModifiedDataNormalized = DataFrame(Dict(Symbol(i) => DFsToConcatNorm[i] for i in 1:NumDataPoints))

    # Remove extreme periods from normalized data before clustering
    NClusters = MinPeriods
    if UseExtremePeriods == 1
        if v println("Pre-removal: ", names(ModifiedDataNormalized)) end
        if v println("Extreme Periods: ", string.(ExtremeWksList)) end
        ClusteringInputDF = select(ModifiedDataNormalized, Not(string.(ExtremeWksList)))
        if v println("Post-removal: ", names(ClusteringInputDF)) end
        NClusters -= length(ExtremeWksList)
    else
        ClusteringInputDF = ModifiedDataNormalized
    end

    ##### STEP 3: Clustering
    cluster_results = []

    # Cluster once regardless of iteration decisions
    push!(cluster_results, cluster(ClusterMethod, ClusteringInputDF, NClusters, nReps, v))

    # Iteratively add worst periods as extreme periods OR increment number of clusters k
    #    until threshold is met or maximum periods are added (If chosen in inputs)
    if (Iterate == 1)
        while (!check_condition(Threshold, last(cluster_results)[1], OldColNames, ScalingMethod, TimestepsPerRepPeriod)) & ((length(ExtremeWksList)+NClusters) < MaxPeriods)
            if IterateMethod == "cluster"
                if v println("Adding a new Cluster! ") end
                NClusters += 1
                push!(cluster_results, cluster(ClusterMethod, ClusteringInputDF, NClusters, nReps, v))
            elseif (IterateMethod == "extreme") & (UseExtremePeriods == 1)
                if v println("Adding a new Extreme Period! ") end
                worst_period_idx = get_worst_period_idx(last(cluster_results)[1])
                removed_period = string(names(ClusteringInputDF)[worst_period_idx])
                select!(ClusteringInputDF, Not(worst_period_idx))
                push!(ExtremeWksList, parse(Int, removed_period))
                if v println(worst_period_idx, " (", removed_period, ") ", ExtremeWksList) end
                push!(cluster_results, cluster(ClusterMethod, ClusteringInputDF, NClusters, nReps, v))
            elseif IterateMethod == "extreme"
                println("INVALID IterateMethod ", IterateMethod, " because UseExtremePeriods is off. Set to 1 if you wish to add extreme periods.")
                break
            else
                println("INVALID IterateMethod ", IterateMethod, ". Choose 'cluster' or 'extreme'.")
                break
            end
        end
        if v && (length(ExtremeWksList)+NClusters == MaxPeriods)
            println("Stopped iterating by hitting the maximum number of periods.")
        elseif v
            println("Stopped by meeting the accuracy threshold.")
        end
    end

    # Interpret Final Clustering Result
    R = last(cluster_results)[1]  # Cluster Object
    A = last(cluster_results)[2]  # Assignments
    W = last(cluster_results)[3]  # Weights
    M = last(cluster_results)[4]  # Centers or Medoids
    DistMatrix = last(cluster_results)[5]  # Pairwise distances
    if v
        println("Total Groups Assigned to Each Cluster: ", W)
        println("Sum Cluster Weights: ", sum(W))
        println("Representative Periods: ", M)
    end

    # K-means/medoids returns indices from DistMatrix as its medoids.
    #   This does not account for missing extreme weeks.
    #   This is corrected retroactively here.
    M = [parse(Int64, string(names(ClusteringInputDF)[i])) for i in M]
    if v println("Fixed M: ", M) end

    ##### Step 4: Aggregation
    # Add the subperiods corresponding to the extreme periods back into the data.
    # Rescale weights to total user-specified number of hours (e.g., 8760 for one year).
    # If LoadExtremePeriod=false (because we don't want to change peak load day), rescale load to ensure total demand is equal

    # Add extreme periods into the clustering result with # of occurences = 1 for each
    ExtremeWksList = sort(ExtremeWksList)
    if UseExtremePeriods == 1
        if v println("Extreme Periods: ", ExtremeWksList) end
        M = [M; ExtremeWksList]
        for w in 1:length(ExtremeWksList)
            insert!(A, ExtremeWksList[w], NClusters+w)
            push!(W, 1)
        end
        NClusters += length(ExtremeWksList) #NClusers from this point forward is the ending number of periods
    end

    N = W  # Keep cluster version of weights stored as N, number of periods represented by RP

    # Rescale weights to total user-specified number of hours
    W = scale_weights(W, WeightTotal, v)

    # Order representative periods chronologically
    #   SORT A W M in conjunction, chronologically by M, before handling them elsewhere to be consistent
    #   A points to an index of M. We need it to point to a new index of sorted M. Hence, AssignMap.
    old_M = M
    df_sort = DataFrame( Weights = W, NumPeriods = N, Rep_Period = M)
    sort!(df_sort, [:Rep_Period])
    W = df_sort[!, :Weights]
    N = df_sort[!, :NumPeriods]
    M = df_sort[!, :Rep_Period]
    AssignMap = Dict( i => findall(x->x==old_M[i], M)[1] for i in 1:length(M))
    A = [AssignMap[a] for a in A]

    # Make PeriodMap, maps each period to its representative period
    PeriodMap = DataFrame(Period_Index = 1:length(A),
                            Rep_Period = [M[a] for a in A],
                            Rep_Period_Index = [a for a in A])

    # Get Symbol-version of column names by type for later analysis
    LoadCols = [Symbol("Load_MW_z"*string(i)) for i in 1:length(load_col_names) ]
    VarCols = [Symbol(var_col_names[i]) for i in 1:length(var_col_names) ]
    FuelCols = [Symbol(fuel_col_names[i]) for i in 1:length(fuel_col_names) ]
    ConstCol_Syms = [Symbol(ConstCols[i]) for i in 1:length(ConstCols) ]

    # Cluster Ouput: The original data at the medoids/centers
    ClusterOutputData = ModifiedData[:,Symbol.(M)]

    # Get zone-wise load multipliers for later scaling in order for weighted-representative-total-zonal load to equal original total-zonal load
    #  (Only if we don't have load-related extreme periods because we don't want to change peak load periods)
    if !LoadExtremePeriod
        load_mults = get_load_multipliers(ClusterOutputData, InputData, M, W, LoadCols, TimestepsPerRepPeriod, NewColNames, NClusters, Ncols)
    end

    # Reorganize Data by Load, Solar, Wind, Fuel, and GrpWeight by Hour, Add Constant Data Back In
    rpDFs = [] # Representative Period DataFrames - Load and Resource Profiles
    gvDFs = [] # Generators Variability DataFrames - Just Resource Profiles
    lpDFs = [] # Load Profile DataFrames - Just Load Profiles
    fpDFs = [] # Fuel Profile DataFrames - Just Fuel Profiles
    for m in 1:NClusters
        rpDF = DataFrame( Dict( NewColNames[i] => ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i] for i in 1:Ncols) )
        gvDF = DataFrame( Dict( NewColNames[i] => ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i] for i in 1:Ncols if (Symbol(NewColNames[i]) in VarCols)) )
        lpDF = DataFrame( Dict( NewColNames[i] => ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i] for i in 1:Ncols if (Symbol(NewColNames[i]) in LoadCols)) )
        if IncludeFuel fpDF = DataFrame( Dict( NewColNames[i] => ClusterOutputData[!,m][TimestepsPerRepPeriod*(i-1)+1 : TimestepsPerRepPeriod*i] for i in 1:Ncols if (Symbol(NewColNames[i]) in FuelCols)) ) end
        if !IncludeFuel fpDF = DataFrame(Placeholder = 1:TimestepsPerRepPeriod) end

        # Add Constant Columns back in
        for c in 1:length(ConstCols)
            rpDF[!,Symbol(ConstCols[c])] .= ConstData[c][1]
            if Symbol(ConstCols[c]) in VarCols
                gvDF[!,Symbol(ConstCols[c])] .= ConstData[c][1]
            elseif Symbol(ConstCols[c]) in FuelCols
                fpDF[!,Symbol(ConstCols[c])] .= ConstData[c][1]
            elseif Symbol(ConstCols[c]) in LoadCols
                lpDF[!,Symbol(ConstCols[c])] .= ConstData[c][1]
            end
        end
        if !IncludeFuel select!(fpDF, Not(:Placeholder)) end

        # Scale Load using previously identified multipliers
        #   Scale lpDF but not rpDF which compares to input data but is not written to file.
        if !LoadExtremePeriod
            for loadcol in LoadCols
                if loadcol ∉ ConstCol_Syms
                    lpDF[!,loadcol] .*= load_mults[loadcol]
                end
            end
        end

        rpDF[!,:GrpWeight] .= W[m]
        rpDF[!,:Cluster] .= M[m]
        push!(rpDFs, rpDF)
        push!(gvDFs, gvDF)
        push!(lpDFs, lpDF)
        push!(fpDFs, fpDF)
    end
    FinalOutputData = vcat(rpDFs...)  # For comparisons with input data to evaluate clustering process
    GVOutputData = vcat(gvDFs...)     # Generators Variability
    LPOutputData = vcat(lpDFs...)     # Load Profiles
    FPOutputData = vcat(fpDFs...)     # Load Profiles

    ##### Step 5: Evaluation

    InputDataTest = InputData[(InputData.Group .<= NumDataPoints*1.0), :]
    ClusterDataTest = vcat([rpDFs[a] for a in A]...) # To compare fairly, load is not scaled here
    RMSE = Dict( c => rmse_score(InputDataTest[:, c], ClusterDataTest[:, c])  for c in OldColNames)

    ##### Step 6: Print to File

    mkpath(joinpath(inpath, TimeDomainReductionFolder))

    ### Load_data_clustered.csv
    load_in = DataFrame(CSV.File(joinpath(inpath, "Load_data.csv"), header=true), copycols=true) #Setting header to false doesn't take the names of the columns; not including it, not including copycols, or, setting copycols to false has no effect
    load_in[!,:Sub_Weights] = load_in[!,:Sub_Weights] * 1.
    load_in[1:length(W),:Sub_Weights] .= W
    load_in[!,:Rep_Periods][1] = length(W)
    load_in[!,:Timesteps_per_Rep_Period][1] = TimestepsPerRepPeriod
    select!(load_in, Not(LoadCols))
    select!(load_in, Not(:Time_Index))
    Time_Index_M = Union{Int64, Missings.Missing}[missing for i in 1:size(load_in,1)]
    Time_Index_M[1:size(LPOutputData,1)] = 1:size(LPOutputData,1)
    load_in[!,:Time_Index] .= Time_Index_M

    for c in LoadCols
        new_col = Union{Float64, Missings.Missing}[missing for i in 1:size(load_in,1)]
        new_col[1:size(LPOutputData,1)] = LPOutputData[!,c]
        load_in[!,c] .= new_col
    end
    load_in = load_in[1:size(LPOutputData,1),:]

    if v println("Writing load file...") end
    CSV.write(joinpath(inpath, Load_Outfile), load_in)

    ### Generators_variability_clustered.csv

    # Reset column ordering, add time index, and solve duplicate column name trouble with CSV.write's header kwarg
    GVColMap = Dict(myinputs["RESOURCE_ZONES"][i] => myinputs["RESOURCES"][i] for i in 1:length(myinputs["RESOURCES"]))
    GVColMap["Time_Index"] = "Time_Index"
    GVOutputData = GVOutputData[!, Symbol.(myinputs["RESOURCE_ZONES"])]
    insertcols!(GVOutputData, 1, :Time_Index => 1:size(GVOutputData,1))
    NewGVColNames = [GVColMap[string(c)] for c in names(GVOutputData)]
    if v println("Writing resource file...") end
    CSV.write(joinpath(inpath, GVar_Outfile), GVOutputData, header=NewGVColNames)

    ### Fuels_data_clustered.csv

    fuel_in = DataFrame(CSV.File(joinpath(inpath, "Fuels_data.csv"), header=true), copycols=true)
    select!(fuel_in, Not(:Time_Index))
    SepFirstRow = DataFrame(fuel_in[1, :])
    NewFuelOutput = vcat(SepFirstRow, FPOutputData)
    rename!(NewFuelOutput, FuelCols)
    insertcols!(NewFuelOutput, 1, :Time_Index => 0:size(NewFuelOutput,1)-1)
    if v println("Writing fuel profiles...") end
    CSV.write(joinpath(inpath, Fuel_Outfile), NewFuelOutput)

    ### Period_map.csv
    if v println("Writing period map...") end
    CSV.write(joinpath(inpath, PMap_Outfile), PeriodMap)

    ### time_domain_reduction_settings.yml
    if v println("Writing .yml settings...") end
    YAML.write_file(joinpath(inpath, YAML_Outfile), myTDRsetup)

    return FinalOutputData, W, RMSE, myTDRsetup, col_to_zone_map
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using CUDArt
using CUDNN
using KUnet
require("mnist.jl")
using MNIST: xtrn, ytrn, xtst, ytst

# load hdf5 file

lr = 0.01

# x1:(28,28,1,64) using 1 channel and 64 minibatch size
@show x1 = Tensor(reshape(MNIST.xtrn[:,1:64], 28, 28, 1, 64))

dims1 = (5,5,1,20)  # using 5x5 filters, 1 input, 20 outputs
conv1 = ConvLayer(w=Filter(float32(randn(dims1)*0.01)),
                  b=Tensor(zeros(Float32,(1,1,dims1[4],1))),
                  pd=PoolingDescriptor((2,2)),
                  pw=UpdateParam(learningRate=lr),
                  pb=UpdateParam(learningRate=2*lr),
                  )

# y1: (24,24,20,64)
# z1: (12,12,20,64) => x2

@show x2=KUnet.forw(conv1,x1)

dims2 = (5,5,20,50)
conv2 = ConvLayer(w=Filter(float32(randn(dims2)*0.01)),
                  b=Tensor(zeros(Float32,(1,1,dims2[4],1))),
                  pd=PoolingDescriptor((2,2)),
                  pw=UpdateParam(learningRate=lr),
                  pb=UpdateParam(learningRate=2*lr),
                  )

# y2: (8,8,50,64)
# z2: (4,4,50,64)
# x3: (800,64)

@show x3=KUnet.forw(conv2,x2)

dims3 = (500,800)
ip3 = Op(w=CudaArray(float32(randn(dims3)*0.01)),
            b=CudaArray(zeros(Float32, dims3[1], 1)),
            f=relu,
            pw=UpdateParam(learningRate=lr),
            pb=UpdateParam(learningRate=2*lr),
            )

# y3: (500,64) => x4

@show x4=KUnet.forw(ip3,x3)

dims4 = (10,500)
ip4 = Op(w=CudaArray(float32(randn(dims4)*0.01)),
            b=CudaArray(zeros(Float32, dims4[1], 1)),
            f=relu,
            pw=UpdateParam(learningRate=lr),
            pb=UpdateParam(learningRate=2*lr),
            )

@show x5=KUnet.forw(ip4,x4)

net4 = [conv1, conv2, ip3, ip4]

ext: jl


text:
# Functions for handling and triggering events

"""
    trigger(event::Event, state, domain=nothing; kwargs...)
    trigger(events::Vector{Event}, state, domain=nothing; kwargs...)

Trigger an `event` or list of `events` if their preconditions hold in the given
`state` and `domain`. See [`execute`](@ref) for keyword arguments.
"""
function trigger(event::Event, state::State,
                 domain::Union{Domain,Nothing}=nothing;
                 as_dist::Bool=false, as_diff::Bool=false)
    # Check whether preconditions hold
    sat, subst = satisfy([event.precond], state, domain; mode=:all)
    if !sat
        @debug "Precondition $(event.precond) does not hold."
        return as_diff ? no_effect(as_dist) : state
    end
    # Update state with effects for each matching substitution
    effects = [substitute(event.effect, s) for s in subst]
    if as_dist
        # Compute product distribution
        eff_dists = [effect_dist(e, state) for e in effects]
        diff = combine(eff_dists...)
    else
        # Accumulate effect diffs
        diff = combine([effect_diff(e, state) for e in effects]...)
    end
    # Return either the difference or the updated state
    return as_diff ? diff : update(state, diff)
end

function trigger(events::Vector{Event}, state::State,
                 domain::Union{Domain,Nothing}=nothing;
                 as_dist::Bool=false, as_diff::Bool=false)
    if length(events) == 0
        diff = as_dist ? DiffDist() : Diff()
    else
        diffs = [trigger(e, state, domain; as_dist=as_dist, as_diff=true)
                 for e in events]
        filter!(d -> !isnothing(d), diffs)
        diff = combine(diffs...)
    end
    # Return either the difference or the updated state
    return as_diff ? diff : update(state, diff)
end


text:
using VoronoiCells
import VoronoiCells: LEFT, RIGHT, LOWER, UPPER
using Base.Test

N = 100
x = 1.0 + rand(N)
y = 1.0 + rand(N)

A = VoronoiCells.voronoiarea(x, y, [LEFT; RIGHT; LOWER; UPPER])
# For most realizations |sum(A) - 1| is approx 10^{-15} and easily
# passing the test. Occasionally it fails with sum(A) = 0.999x, where x != 9
@test_approx_eq sum(A) 1.0

# ------------------------------------------------------------
# Compare with Deldir

using Deldir
A2 = Deldir.voronoiarea(x, y; rw=[LEFT; RIGHT; LOWER; UPPER])
@test_approx_eq A A2


text:
__precompile__()
module eps_InAs_ptype
using QuadGK
export epsFCL_ptype,Gamma_ptype,m_star_ptype
#from integration import Integrate
#from scipy.integrate import quad
 #but I can't use QuadGk for larger values of eV, can only functionine function Integrate, tested E = 0.0037eV, get same eps of 10.240628772514803 + 0.03160036047883907im

const e = 1.602*^(10,-19)
const eV = 1.60218*^(10,-19)
const eps_0 = 8.8542*^(10,-12)
const c = 3*^(10,8)
const hbar = 1.05457*^(10,-34) #m2kg/s
const hbar_eV = 6.5821*^(10,-16) #eV s
const m_e =9.109*^(10,-31)#kg
const m_0 = 0.9109*^(10,-30)
const m_0_eV = 0.51099906*^(10,6) #eV

const E0_InAs = 0.42
const dt0_InAs = 0.28*0.001
const b0_InAs = 93
const A_InAs = 14*10^(-9) #for n-type  InAs
const B_InAs = 1.97*10^(-7)
const C_InAs = 57.9*10^(-12) #N in m^(-3)

const P_square = 11.9   #eV
const E_g_InAs= 0.416 #eV

const umin =20.0
const umax = 530.0
const Nref = 1.1*10^17*10^6 #in m^(-3) since N_base in m^(-3)
const phi = 0.46
const t1 = 2.3
const t2 = 3.0

@inline function m_star_ptype(n) # change m_e according to MB shift
    term = 8.0*P_square/eV*hbar^2*(3.0*pi^2*n)^(2/3)/(3.0*m_e*E_g_InAs^2)
    return (1.0+(4.0*P_square/(3.0*E_g_InAs))*(1.0+term)^(-1/2))^(-1)*m_0
end
precompile(m_star_ptype, (Float64,))

@inline function Gamma_ptype(N_Base,T)
    mu_franc =  umin + (umax*(300.0/T)^t1-umin)/(1+(N_Base/(Nref*(T/300.0)^t2))^phi)#from Francoeur
    return e/(m_star_ptype(N_Base)*mu_franc*10000.0) #s^(-1) change mu from cm^(-2) to m^(-2)
end
precompile(Gamma_ptype, (Float64,Float64))

const eps_inf_InAs = 12.3
const g = 2.89*10.0^11  #s^(-1), gamma
const o_TO = 4.12*10.0^(13) #s^(-1)
const o_LO = 4.58*10.0^(13)
@inline function epsFCL_ptype(omega,mstar,gamma,N_base)
    o_p_square=N_base*e^2/(12.3*eps_0*mstar) #s^(-2)
    return eps_inf_InAs*(1.0+(o_LO^2-o_TO^2)/(o_TO^2-omega^2-im*omega*g)-o_p_square/(omega*(omega+ im*gamma)))
end
precompile(epsFCL_ptype,(Float64,Float64,Float64,Float64))
end
#@inline function epsIB(omega,N0,T,E0_T_InAs_value) #use the function from ntype
#    if(hbar_eV*omega<=E0_T_InAs_value)
#        alpha_lb_2=0.0
#    else
#        alpha_lb_2= (3*10^8*(hbar_eV*omega-E0_T_InAs_value))^(1/2)*100
#    end
#    m_IB_pp = alpha_lb_2*c/(2*omega)

#    integratedfn(omega_p) = alpha_lb_2(omega_p,N0,T)/(omega^2-omega_p^2)
#    if(hbar_eV*omega<E0_T_InAs_value)
#        m_IB_p_value =1.0
#    else
     #   result  =1.0 + (c/pi)*quadgk(integratedfn,100.0,10.0^25,rtol=1e-8)[1] # result  =(1.0 + (c/pi)*Integrate(integratedfn, 10000,100.0,10.0^25))
#        m_IB_p_value = (1.0 + 0.00026172525)
#    end
#    return (m_IB_p_value^2-m_IB_pp^2) + im*2*m_IB_p_value*m_IB_pp
#end
#precompile(epsIB,(Float64,Float64,Float64,Float64))

#function eps_InAs_ntype(E_photon::Float64,N_base::Float64)::ComplexF64 #N_base in m^(-3)
    #E_photon= eVtoOmega(omega)
#    omega = eVtoOmega(E_photon)
#    T=300.0 #o_angular = 5.38*10^14   #hbar*5.38*10^14/eV = 0.35eV
#    if(E_photon>0.35)
#        return epsIB(omega,N_base,T) +  epsFCL(omega,N_base,T)
#    else
#        return epsFCL(omega,N_base,T)
#    end
#end

#omegalist= 2.0*10.0^14:5.75*10.0^13:2.5*10.0^15#
#eVlist=  0.2:0.02:1.5
#epslist = []
#for ind = 1:length(eVlist)
#    eps = eps_InAs_ntype(eVlist[ind],10.0^23)
#    println("eps_InAs",eps)
#    push!(epslist, imag(eps))
#end
#plot(eVlist,epslist)
#savefig("test_InAs_real.png")


text:
using Koala
using Revise
using KoalaTrees
using Base.Test

# test the histogram maker:
@test KoalaTrees.histogram(21:41, 11:1000, 3)[1] ==
[[21, 22, 23],
 [24, 25],
 [26, 27, 28],
 [29, 30],
 [31, 32, 33],
 [34, 35],
 [36, 37, 38],
 [39, 40, 41]]

# test clean
TreeRegressor(regularization=1.0)
TreeRegressor(min_patterns_split=1)
TreeRegressor(max_bin=1)

const X, y = load_ames();

const all = eachindex(y) # iterator for all rows

const train, test = split(all, 0.8); # 80:20 split

const rgs = TreeRegressor(penalty=0.5)

const mach = Machine(rgs, X, y, train, features=names(X)[1:end-1])
fit!(mach, train)
showall(mach)
score = err(mach, test)
println("error = $score")
@test score > 4e4 && score < 5e4
score = mean(cv(mach, all))
println("error = $score")

rgs.extreme = true
fit!(mach)
@test err(mach, test) < 5e4

rgs.penalty = 0
rgs.extreme = false
rgs.max_bin = 16
rgs.bin_factor = 1
fit!(mach)
err(mach, test)
@test err(mach, test) < 5e4


text:
using GenCurveSegmentation
using Test

@testset "GenCurveSegmentation.jl" begin
    # Write your tests here.
end


text:
export
    geodesic_distance,
    geographical_to_geocentric_lat

"""
    geodesic_distance(L1, B1, L2, B2)

Computes the geodesic distance between two points on the Earth with a ~50 meters precision

Parameters:
L1, B1 : longitude and latitude of the first point (in radians)
L2, B2 : longitude and latitude of the second point (in radians)

Returns:
    distance in km
"""
function geodesic_distance(L1, B1, L2, B2)
    global earth_equ_radius, earth_flattening
    F = (B1+B2)/2
    G = (B1-B2)/2
    lambda = (L1-L2)/2
    sF = sin(F)
    cF = cos(F)
    sG = sin(G)
    cG = cos(G)
    sl = sin(lambda)
    cl = cos(lambda)

    S = sG*sG*cl*cl+cF*cF*sl*sl
    C = cG*cG*cl*cl+sF*sF*sl*sl

    omega = atan(sqrt(S/C))
    R = sqrt(S*C)/omega
    H1 = (3*R - 1)/2/C
    H2 = (3*R + 1)/2/S
    D  = 2*omega*earth_equ_radius

    return D * (1 + earth_flattening * (H1 * sF * sF * cG * cG - H2 * cF * cF * sG *sG))
end

"""
    geographical_to_geocentric_lat(phi, H)

Convert geographical latitude to geocentric latitude

[Meeus-1998: chapter 11]

Parameters:
    phi : geographical latitude in radians
    H   : altitude of the observer above sea level

Returns:
    phi1: geocentric latitude in radians
    rho*sin(phi1)
    rho*cos(phi1)
"""
function geographical_to_geocentric_lat(phi, H)
    global earth_flattening
    ratio = (1-earth_flattening)
    u = atan(ratio*sin(phi), cos(phi))
    r = H/6378140
    rhosinphi1 = ratio*sin(u) + r*sin(phi)
    rhocosphi1 = cos(u) + r*cos(phi)
    phi1 = atan(rhosinphi1, rhocosphi1)
    return (phi1, rhosinphi1, rhocosphi1)
end


text:
# baremodule Defaults
# using CoilFunctions, Unitful
# # using Unitful
#
# # abstract type Product end
# # struct SBE4 <: Product end
# # struct SBE5 <: Product end
# # struct SBE7 <: Product end
# # struct SBE9 <: Product end
# #
# # CoilGeometry(::Type{SBE4}) = CoilGeometry(66.5u"mm",83.5u"mm",21u"mm")
# # CoilGeometry(::Type{SBE5}) = CoilGeometry(86.5u"mm",106u"mm",26u"mm")
# # CoilGeometry(::Type{SBE7}) = CoilGeometry(104u"mm",136u"mm",26u"mm")
# # CoilGeometry(::Type{SBE9}) = CoilGeometry(144u"mm",186u"mm",26u"mm")
# # CoilGeometry(::SBE5) = CoilGeometry(86.5u"mm",106u"mm",26u"mm")
# #
# # power(::Type{SBE4})=30.0u"W"
# # power(::Type{SBE5})=44.0u"W"
# # power(::Type{SBE7})=50.0u"W"
# # power(::Type{SBE9})=86.0u"W"
#
# # coil=Dict(:SBE4 => CoilGeometry(66.5u"mm",83.5u"mm",21u"mm"),
# #     :SBE5 => CoilGeometry(86.5u"mm",106u"mm",26u"mm"),
# #     :SBE7 => CoilGeometry(104u"mm",136u"mm",26u"mm"),
# #     :SBE9 => CoilGeometry(144u"mm",186u"mm",26u"mm"))
#
# end

const coil=Dict(:SBE4 => CoilGeometry(66.5u"mm",83.5u"mm",21u"mm"),
    :SBE5 => CoilGeometry(86.5u"mm",106u"mm",26u"mm"),
    :SBE7 => CoilGeometry(104.0u"mm",136u"mm",26u"mm"),
    :SBE9 => CoilGeometry(144.0u"mm",186u"mm",26u"mm"))

const coilcavity=Dict(:SBE4 => CoilGeometry(62.5u"mm",83.5u"mm",25u"mm"),
    :SBE5 => CoilGeometry(82.5u"mm",110u"mm",30u"mm"),
    :SBE7 => CoilGeometry(100.0u"mm",140u"mm",30u"mm"),
    :SBE9 => CoilGeometry(140.0u"mm",190u"mm",30u"mm"))

const power=Dict(:SBE4 => 30.0u"W", :SBE5 => 44.0u"W", :SBE7=>50.0u"W",:SBE9=>86.0u"W")
Ω(s::Symbol,v::Unitful.Voltage) = v^2/power[s] |> u"Ω"


text:
export RQ, rq, rq!

"""
    RQ <: Factorization

An RQ matrix factorization stored in a packed format, typically obtained from
[`rq`](@ref). If ``A`` is an `m`×`n` matrix, then

```math
A = [0 R] Q
```

where ``Q`` is an orthogonal/unitary matrix and ``R`` is upper triangular (trapezoidal
if `m < n`).
The matrix ``Q`` is stored as a sequence of Householder reflectors ``v_i``
and coefficients ``\\tau_i`` where:

```math
Q = \\prod_{i=1}^{\\min(m,n)} (I - \\tau_i v_i v_i^T).
```

Iterating the decomposition produces the components `R` and `Q`.

The object has two fields:

* `factors` is an `m`×`n` matrix.

  - The upper triangular part contains the elements of ``R``, that is `R =
    triu(F.factors)` for a `QR` object `F`.

  - The subdiagonal part contains the reflectors ``v_i`` stored in a packed format where
    ``v_i`` is the ``i``th column of the matrix `V = I + tril(F.factors, -1)`.

* `τ` is a vector  of length `min(m,n)` containing the coefficients ``\tau_i``.

"""
struct RQ{T,S<:AbstractMatrix{T}} <: Factorization{T}
    factors::S
    τ::Vector{T}

    function RQ{T,S}(factors, τ) where {T,S<:AbstractMatrix{T}}
        require_one_based_indexing(factors)
        new{T,S}(factors, τ)
    end
end
RQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = RQ{T,typeof(factors)}(factors, τ)
function RQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T}
    RQ(convert(AbstractMatrix{T}, factors), convert(Vector{T}, τ))
end

# iteration for destructuring into components
Base.iterate(S::RQ) = (S.R, Val(:Q))
Base.iterate(S::RQ, ::Val{:Q}) = (S.Q, Val(:done))
Base.iterate(S::RQ, ::Val{:done}) = nothing

Base.size(F::RQ, dim::Integer) = size(getfield(F, :factors), dim)
Base.size(F::RQ) = size(getfield(F, :factors))

function Base.getproperty(F::RQ, d::Symbol)
    m, n = size(F)
    if d === :R
        if m <= n
            return triu!(getfield(F, :factors)[1:m, n-m+1:n])
        else
            return triu(getfield(F, :factors),min(0,n-m))
        end
    elseif d === :Q
        return RQPackedQ(getfield(F, :factors), F.τ)
    else
        getfield(F, d)
    end
end

"""
    RQPackedQ <: AbstractMatrix

The orthogonal/unitary ``Q`` matrix of an RQ factorization stored in [`RQ`](@ref) format.
"""
struct RQPackedQ{T,S<:AbstractMatrix{T}} <: AbstractQ{T}
    factors::S
    τ::Vector{T}

    function RQPackedQ{T,S}(factors, τ) where {T,S<:AbstractMatrix{T}}
        require_one_based_indexing(factors)
        new{T,S}(factors, τ)
    end
end

# Conversions
RQ{T}(A::RQ) where {T} = RQ(convert(AbstractMatrix{T}, A.factors), convert(AbstractVector{T}, A.τ))
Factorization{T}(A::RQ{T}) where {T} = A
Factorization{T}(A::RQ) where {T} = RQ{T}(A)
AbstractMatrix(F::RQ) = F.R * F.Q
AbstractArray(F::RQ) = AbstractMatrix(F)
Matrix(F::RQ) = Array(AbstractArray(F))
Array(F::RQ) = Matrix(F)

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, F::RQ)
    summary(io, F); println(io)
    println(io, "R factor:")
    show(io, mime, F.R)
    println(io, "\nQ factor:")
    show(io, mime, F.Q)
end

RQPackedQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = RQPackedQ{T,typeof(factors)}(factors, τ)
function RQPackedQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T}
    RQPackedQ(convert(AbstractMatrix{T}, factors), convert(Vector{T}, τ))
end
RQPackedQ{T}(Q::RQPackedQ) where {T} = RQPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ))
AbstractMatrix{T}(Q::RQPackedQ{T}) where {T} = Q
AbstractMatrix{T}(Q::RQPackedQ) where {T} = RQPackedQ{T}(Q)
# this makes the rectangular version, alas.
# function Matrix{T}(Q::RQPackedQ{T}) where {T<:BlasFloat}
#     m,n = size(Q.factors)
#     i0= m > n ? (m-n+1) : 1
#     LAPACK.orgrq!(Q.factors[i0:m,1:n], Q.τ)
# end

Matrix{T}(Q::RQPackedQ{S}) where {T,S} = Matrix{T}(lmul!(Q, Matrix{S}(I, size(Q,2), size(Q,2))))
Base.size(Q::RQPackedQ, dim::Integer) = size(getfield(Q, :factors), dim == 1 ? 2 : dim)

function lmul!(A::RQPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix}
    m,n = size(A.factors)
    if m > n
        return LAPACK.ormrq!('L','N',view(A.factors, m-n+1:m, 1:n),A.τ,B)
    else
        return LAPACK.ormrq!('L','N',A.factors,A.τ,B)
    end
end
function lmul!(A::RQPackedQ, B::AbstractVecOrMat)
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if nA != mB
        throw(DimensionMismatch("matrix A has dimensions ($nA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    mm = min(mA,nA)
    @inbounds begin
        for k = mm:-1:1
            krow = mA-mm+k
            kcol = nA-mm+k
            for j = 1:nB
                vBj = B[kcol,j]
                for i = 1:kcol-1
                    # vBj += conj(Afactors[i,k])*B[i,j]
                    vBj += (Afactors[krow,i])*B[i,j]
                end
                vBj = conj(A.τ[k])*vBj
                B[kcol,j] -= vBj
                for i = 1:kcol-1
                    # B[i,j] -= Afactors[i,k]*vBj
                    B[i,j] -= conj(Afactors[krow,i])*vBj
                end
            end
        end
    end
    B
end

function (*)(B::StridedMatrix, Q::RQPackedQ)
    TBQ = promote_type(eltype(Q), eltype(B))
    Qnew = convert(AbstractMatrix{TBQ}, Q)
    if size(Q.factors, 2) == size(B, 2)
        Bnew = copy_oftype(B, TBQ)
    elseif size(Q.factors, 1) == size(B, 2)
        # People will expect this, misguided though it seems.
        Bnew = hcat(zeros(TBQ, size(B,1), size(Q.factors, 2) - size(B,2)),B)
    else
        throw(DimensionMismatch("first dimension of matrix must have size either $(size(Q.factors, 1)) or $(size(Q.factors, 2))"))
    end
    rmul!(Bnew, Qnew)
end

function lmul!(adjA::Adjoint{<:Any,<:RQPackedQ{T,S}}, B::StridedVecOrMat{T}
               ) where {T<:BlasReal,S<:StridedMatrix}
    A = adjA.parent
    m,n = size(A.factors)
    if m > n
        return LAPACK.ormrq!('L','T',view(A.factors, m-n+1:m, 1:n),A.τ,B)
    else
        return LAPACK.ormrq!('L','T',A.factors,A.τ,B)
    end
end
function lmul!(adjA::Adjoint{<:Any,<:RQPackedQ{T,S}}, B::StridedVecOrMat{T}
               ) where {T<:BlasComplex,S<:StridedMatrix}
    A = adjA.parent
    m,n = size(A.factors)
    if m > n
        return LAPACK.ormrq!('L','C',view(A.factors, m-n+1:m, 1:n),A.τ,B)
    else
        return LAPACK.ormrq!('L','C',A.factors,A.τ,B)
    end
end
function lmul!(adjA::Adjoint{<:Any,<:RQPackedQ}, B::AbstractVecOrMat)
    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB, nB = size(B,1), size(B,2)
    if nA != mB
        throw(DimensionMismatch("matrix A has dimensions ($nA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors
    mm = min(mA,nA)
    @inbounds begin
        for k = 1:mm
            krow = mA-mm+k
            kcol = nA-mm+k
            for j = 1:nB
                vBj = B[kcol,j]
                for i = 1:kcol-1
                    # vBj += conj(Afactors[i,k])*B[i,j]
                    vBj += (Afactors[krow,i])*B[i,j]
                end
                vBj = A.τ[k]*vBj
                B[kcol,j] -= vBj
                for i = 1:kcol-1
                    # B[i,j] -= Afactors[i,k]*vBj
                    B[i,j] -= conj(Afactors)[krow,i]*vBj
                end
            end
        end
    end
    B
end
function rmul!(A::StridedVecOrMat{T}, B::RQPackedQ{T,S}
               ) where {T<:BlasFloat,S<:StridedMatrix}
    m,n = size(B.factors)
    if m > n
        return LAPACK.ormrq!('R','N',view(B.factors, m-n+1:m, 1:n),B.τ,A)
    else
        return LAPACK.ormrq!('R','N',B.factors,B.τ,A)
    end
end
function rmul!(A::StridedMatrix,Q::RQPackedQ)
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != nQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($nQ, $nQ)"))
    end
    Qfactors = Q.factors
    mm = min(mQ, nQ)
    @inbounds begin
        for k = 1:min(mQ,nQ)
            krow = mQ-mm+k
            kcol = nQ-mm+k
            for i = 1:mA
                vAi = A[i,kcol]
                for j = 1:kcol-1
                    # vAi += A[i,j]*Qfactors[j,k]
                    vAi += A[i,j]*conj(Qfactors[krow,j])
                end
                vAi = vAi*conj(Q.τ[k])
                A[i,kcol] -= vAi
                for j = 1:kcol-1
                    # A[i,j] -= vAi*conj(Qfactors[j,k])
                    A[i,j] -= vAi*(Qfactors[krow,j])
                end
            end
        end
    end
    A
end
function rmul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:RQPackedQ{T}}
               ) where {T<:BlasReal}
    B = adjB.parent
    m,n = size(B.factors)
    if m > n
        return LAPACK.ormrq!('R','T',view(B.factors, m-n+1:m, 1:n),B.τ,A)
    else
        return LAPACK.ormrq!('R','T',B.factors,B.τ,A)
    end
end
function rmul!(A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:RQPackedQ{T}}
               ) where {T<:BlasComplex}
    B = adjB.parent
    m,n = size(B.factors)
    if m > n
        return LAPACK.ormrq!('R','C',view(B.factors, m-n+1:m, 1:n),B.τ,A)
    else
        return LAPACK.ormrq!('R','C',B.factors,B.τ,A)
    end
end
function rmul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:RQPackedQ})
    Q = adjQ.parent
    mQ, nQ = size(Q.factors)
    mA, nA = size(A,1), size(A,2)
    if nA != nQ
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but matrix Q has dimensions ($nQ, $nQ)"))
    end
    Qfactors = Q.factors
    mm = min(mQ, nQ)
    @inbounds begin
        for k = min(mQ,nQ):-1:1
            krow = mQ-mm+k
            kcol = nQ-mm+k
            for i = 1:mA
                vAi = A[i,kcol]
                for j = 1:kcol-1
                    # vAi += A[i,j]*Qfactors[j,k]
                    vAi += A[i,j]*conj(Qfactors[krow,j])
                end
                vAi = vAi*Q.τ[k]
                A[i,kcol] -= vAi
                for j = 1:kcol-1
                    # A[i,j] -= vAi*conj(Qfactors[j,k])
                    A[i,j] -= vAi*(Qfactors[krow,j])
                end
            end
        end
    end
    A
end

ldiv!(A::RQ{T}, B::StridedMatrix{S}) where {T,S} = lmul!(adjoint(A.Q), ldiv!(UpperTriangular(A.R), B))
ldiv!(A::RQ{T}, B::StridedVector{S}) where {T,S} = lmul!(adjoint(A.Q), ldiv!(UpperTriangular(A.R), B))

# apply shifted reflector from right
@inline function reflectorYlppa!(x::AbstractVector, τ::Number, A::StridedMatrix)
    require_one_based_indexing(x)
    m, n = size(A)
    if length(x) != n
        throw(DimensionMismatch("reflector has length $(length(x)), which must match the second dimension of matrix A, $n"))
    end
    @inbounds begin
        for j = 1:m
            # dot
            Ajv = A[j, n]
            for i = 1:n-1
                Ajv += A[j, i] * x[i+1]
            end

            Ajv = τ * Ajv

            # ger
            A[j, n] -= Ajv
            for i = 1:n-1
                A[j, i] -= Ajv * x[i+1]'
            end
        end
    end
    return A
end

function _rqfactUnblocked!(A::AbstractMatrix{T}) where {T}
    require_one_based_indexing(A)
    m, n = size(A)
    τ = zeros(T, min(m,n))
    v = zeros(T, max(m,n))
    mm = min(m,n)
    for k = mm:-1:1
        krow = m-mm+k
        kcol = n-mm+k
        v[1] = A[krow,kcol]'
        for i = 2:kcol
            v[i] = A[krow,i-1]'
        end
        x = view(v, 1:kcol)
        if kcol==1 && T<:Real
            τk = zero(T)
        else
            τk = reflector!(x)
        end
        τ[k] = τk
        reflectorYlppa!(x, τk, view(A, 1:krow-1, 1:kcol))
        A[krow,kcol] = v[1]'
        for i = 2:kcol
            A[krow,i-1] = v[i]'
        end
    end
    RQ(A, τ)
end

"""
    rq(A) -> F

Compute the RQ factorization of the matrix `A`: an orthogonal/unitary matrix `Q` and
and upper triangular or trapezoidal matrix `R` such that
```math
A = [0 R] Q
```
The returned object `F` stores the factorization in a packed format [`RQ`](@ref).
Unlike the `QR` factorization this does not handle least squares problems.
If `A` is `m`×`n` then `Q` is `n`×`n` and `R` is `m`×`min(m,n)`.
"""
function rq end

"""
    rq!(A) -> F

The same as [`rq`](@ref), but saves space by overwriting the input `A`.
"""
rq!(A::StridedMatrix{<:BlasFloat}) = RQ(LAPACK.gerqf!(A)...)
rq!(A::StridedMatrix) = _rqfactUnblocked!(A)

rq(x::Number) = rq(fill(x,1,1))
function rq(v::AbstractVector)
    require_one_based_indexing(v)
    rq(reshape(v, (length(v), 1)))
end

# There is a local version of this:
# using LinearAlgebra: _qreltype

function rq(A::AbstractMatrix{T}) where {T}
    require_one_based_indexing(A)
    AA = similar(A, LinearAlgebra._qreltype(T), size(A))
    copyto!(AA, A)
    return rq!(AA)
end


text:
using DarkTracks
using Test

videofile = joinpath(@__DIR__(), "..", "test", "testvideo.mp4")
ts, inds = main(videofile, stop_frame = 100, spatialscale = 10)

@test ts == [1.32, 1.54, 1.76, 1.98]
@test inds == [(510, 730), (520, 730), (520, 730), (520, 730)]


text:
"""
# AbstractThing
The `AbstractThing` struct is placed at the top of the taxonomy.
It makes it possible to define pseudo-universal functions to facilitate the manipulation of objects in general.
"""
abstract type AbstractThing end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
@everywhere using Roadnet_MDP
@everywhere importall POMDPs, POMDPToolbox
@everywhere using MCTS
@everywhere using DiscreteValueIteration
@everywhere include("roadnet_pursuer_MDP.jl")
# for parallel policy search, use julia -p 2 (for two threads)
using JLD
using PyPlot

include("network_library.jl")

ext_rwd = 2000.
cgt_rwd = -2000.
sen_rwd = -100.
go = original_roadnet(exit_rwd=ext_rwd,caught_rwd=cgt_rwd,sensor_rwd=sen_rwd)
gm = medium_roadnet(exit_rwd=ext_rwd,caught_rwd=cgt_rwd,sensor_rwd=sen_rwd)

fold = "mcts_vi_comparison_data/" #folder for data/images to go

display_network(go,evader_locs=[1],pursuer_locs=[4],scale=1.0,fname="$(fold)go_net",ftype=:svg)
display_network(gm,evader_locs=[1],pursuer_locs=[4],scale=1.0,fname="$(fold)gm_net",ftype=:svg)

t = 0.80
d = 0.95

mdp_go = roadnet_with_pursuer(go,tp=t,d=d)
mdp_gm = roadnet_with_pursuer(gm,tp=t,d=d)

pygui(false)

function plot_results()
    results_set = jldopen("$(fold)mcts_vi_comparison_results_set.jld","r") do file
	read(file,"results_set")
	end

    fig, ax = PyPlot.subplots(2,1)
    i = 1
    for (net,res) in results_set
        #  display(res)
        ax[i][:violinplot](res,widths=0.5,points=500,showmedians=true)
        i+=1
    end
    PyPlot.savefig("$(fold)mcts_vi_comparison_plot.pdf",transparent=true)
     # the :gm plot isn't working very well, and I don't think I need it, I can just use the :go graph as an example

    figgo,axgo = PyPlot.subplots(1,1)
    figgo[:set_size_inches](10,6)
    fsize = 15

    solve_set = [1,2,3,4,7]
    solve_labels = ["MCTS\nd=1","MCTS\nd=3","MCTS\nd=6","MCTS\nd=10","MCTS\nd=15","MCTS\nd=30","Value\nIteration"]
    s_ticks = collect(1:length(solve_set))
    axgo[:violinplot](results_set[:go][:,solve_set],widths=0.5,points=500,showmedians=true)
    axgo[:set_xlabel]("Solver",fontsize=fsize)
    axgo[:set_ylabel]("Expected Rewards",fontsize=fsize)
    axgo[:set_title]("Comparison of MCTS policies to Value Iteration policy",fontsize=fsize)
    axgo[:set_xticks](s_ticks)
    axgo[:xaxis][:set_ticklabels](solve_labels[solve_set])
    fig[:tight_layout]()
    PyPlot.savefig("$(fold)mcts_vi_comparison_plot_orig_net.pdf",transparent=true)
end

function get_results(mdp_go,mdp_gm;num_repeats=5,display_matrix=false)
    vi_policies = jldopen("$(fold)mcts_vi_comparison_valueiteration.jld","r") do file
	read(file,"solns")
	end
    vi_solvers = jldopen("$(fold)mcts_vi_comparison_valueiteration.jld","r") do file
	read(file,"solvers")
	end

    mcts_policies = jldopen("$(fold)mcts_vi_comparison_mcts.jld","r") do file
	read(file,"solns")
	end
    mcts_solvers = jldopen("$(fold)mcts_vi_comparison_mcts.jld","r") do file
	read(file,"solvers")
	end

    policy_set = Dict(:vi=>vi_policies,:mcts=>mcts_policies)
    sovler_set = Dict(:vi=>vi_solvers,:mcts=>mcts_solvers)
    mdp_set = Dict(:go=>mdp_go,:gm=>mdp_gm)

    num_vi_policies = length(vi_policies[:go])
    num_mcts_policies = length(mcts_policies[:go])
    num_po = num_vi_policies+num_mcts_policies
    results_set = Dict(:go=>zeros(num_repeats,num_po),:gm=>zeros(num_repeats,num_po))

    for (solver_type,mdp_policies) in policy_set
        println(solver_type)
        for (mdp_type,policies) in mdp_policies
            println(mdp_type)
            for (pnum,p) in enumerate(policies)
                println("Running policy $pnum of $solver_type/$mdp_type")
                # simulate performance of policies on the appropraite mdp
                max_steps = 100

                PT = repmat([p],1,num_repeats)
                initial_state = roadnet_pursuer_state(1,4)

                q = []
                for pol in PT
                    value = Sim(mdp_go,pol,initial_state,max_steps=max_steps)
                    push!(q,value)
                end
                sim_results = run_parallel(q) do sim,hist
                    #  display(discounted_reward(hist))
                    return [:steps=>n_steps(hist), :reward=>undiscounted_reward(hist)]
                end

                rewards = [float(x) for x in sim_results[:reward]] #need to convert to float, because sim_results[:reward] is a union of floats and missing values, even though there aren't any missing values...
                if solver_type == :mcts
                    results_set[mdp_type][:,pnum] = rewards
                else
                    results_set[mdp_type][:,end] = rewards
                end

                if display_matrix
                    println(size(results_set[mdp_type]))
                    display(results_set[mdp_type])
                end
            end
        end
    end

    jldopen("$(fold)mcts_vi_comparison_results_set.jld","w") do file
        write(file,"results_set",results_set)
    end

    return results_set
end

function get_policies(mdp_go,mdp_gm;solve_vi=false)
	its = 2000
	exp_const = 2000.0

	println("Initializing Sovlers")

	# mcts
	s_mcts_30 = MCTSSolver(n_iterations=its,depth=30,exploration_constant=exp_const,enable_tree_vis=true)
	s_mcts_15 = MCTSSolver(n_iterations=its,depth=15,exploration_constant=exp_const,enable_tree_vis=true)
	s_mcts_10 = MCTSSolver(n_iterations=its,depth=10,exploration_constant=exp_const,enable_tree_vis=true)
	s_mcts_6 = MCTSSolver(n_iterations=its,depth=6,exploration_constant=exp_const,enable_tree_vis=true)
	s_mcts_3 = MCTSSolver(n_iterations=its,depth=3,exploration_constant=exp_const,enable_tree_vis=true)
	s_mcts_1 = MCTSSolver(n_iterations=its,depth=1,exploration_constant=exp_const,enable_tree_vis=true)
    mcts_solvers = Dict(:mcts=>[s_mcts_10,s_mcts_6,s_mcts_3,s_mcts_1])

	# value iteration
	s_vi = ValueIterationSolver(max_iterations=150,belres=1e-6)
    vi_solvers = Dict(:vi=>[s_vi])

	if solve_vi
		# Value iteration takes FOREVER. Don't do this everytime
		println("Running ValueIterationSolver on original roadnet")
		vi_orig_soln = solve(s_vi,mdp_go)
		println("Running ValueIterationSolver on medium roadnet")
		vi_med_soln = solve(s_vi,mdp_gm)

        solns = Dict(:gm=>[vi_med_soln],:go=>[vi_orig_soln])

        jldopen("$(fold)mcts_vi_comparison_valueiteration.jld","w") do file
			addrequire(file,DiscreteValueIteration)
			write(file,"solns",solns)
            write(file,"solvers",vi_solvers)
		end
	end

	println("Running MCTS")
	mcts_orig_30_soln = solve(s_mcts_30,mdp_go)
	mcts_med_30_soln = solve(s_mcts_30,mdp_gm)

	mcts_orig_15_soln = solve(s_mcts_15,mdp_go)
	mcts_med_15_soln = solve(s_mcts_15,mdp_gm)

	mcts_orig_10_soln = solve(s_mcts_10,mdp_go)
	mcts_med_10_soln = solve(s_mcts_10,mdp_gm)

	mcts_orig_6_soln = solve(s_mcts_6,mdp_go)
	mcts_med_6_soln = solve(s_mcts_6,mdp_gm)

	mcts_orig_3_soln = solve(s_mcts_3,mdp_go)
	mcts_med_3_soln = solve(s_mcts_3,mdp_gm)

	mcts_orig_1_soln = solve(s_mcts_1,mdp_go)
	mcts_med_1_soln = solve(s_mcts_1,mdp_gm)

    mcts_solns = Dict(:gm=>[mcts_med_1_soln,mcts_med_3_soln,mcts_med_6_soln,mcts_med_10_soln,mcts_med_15_soln,mcts_med_30_soln],:go=>[mcts_orig_1_soln,mcts_orig_3_soln,mcts_orig_6_soln,mcts_orig_10_soln,mcts_orig_15_soln,mcts_orig_30_soln])

	# write the policies to a file
    jldopen("$(fold)mcts_vi_comparison_mcts.jld","w") do file
		addrequire(file,DiscreteValueIteration)
		write(file,"solns",mcts_solns)
        write(file,"solvers",mcts_solvers)
	end
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
# finite mixture models

"""

  All subtypes of `AbstractMixtureModel` should implement the following methods:

  - ncomponents(d): the number of components

  - component(d, k):  return the k-th component

  - probs(d):       return a vector of prior probabilities over components.
"""
abstract type AbstractMixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution} <: Distribution{VF, VS} end

"""
MixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution,CT<:Real}
A mixture of distributions, parametrized on:
* `VF,VS` variate and support
* `C` distribution family of the mixture
* `CT` the type for probabilities of the prior
"""
struct MixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution,CT<:Real} <: AbstractMixtureModel{VF,VS,C}
    components::Vector{C}
    prior::Categorical{CT}

    function MixtureModel{VF,VS,C}(cs::Vector{C}, pri::Categorical{CT}) where {VF,VS,C,CT}
        length(cs) == ncategories(pri) ||
            error("The number of components does not match the length of prior.")
        new{VF,VS,C,CT}(cs, pri)
    end
end

const UnivariateMixture{S<:ValueSupport,   C<:Distribution} = AbstractMixtureModel{Univariate,S,C}
const MultivariateMixture{S<:ValueSupport, C<:Distribution} = AbstractMixtureModel{Multivariate,S,C}
const MatrixvariateMixture{S<:ValueSupport,C<:Distribution} = AbstractMixtureModel{Matrixvariate,S,C}

# Interface

"""
    component_type(d::AbstractMixtureModel)

The type of the components of `d`.
"""
component_type(d::AbstractMixtureModel{VF,VS,C}) where {VF,VS,C} = C

"""
    components(d::AbstractMixtureModel)

Get a list of components of the mixture model `d`.
"""
components(d::AbstractMixtureModel)

"""
    probs(d::AbstractMixtureModel)

Get the vector of prior probabilities of all components of `d`.
"""
probs(d::AbstractMixtureModel)

"""
    mean(d::Union{UnivariateMixture, MultivariateMixture})

Compute the overall mean (expectation).
"""
mean(d::AbstractMixtureModel)

"""
    insupport(d::MultivariateMixture, x)

Evaluate whether `x` is within the support of mixture distribution `d`.
"""
insupport(d::AbstractMixtureModel, x::AbstractVector)

"""
    pdf(d::Union{UnivariateMixture, MultivariateMixture}, x)

Evaluate the (mixed) probability density function over `x`. Here, `x` can be a single
sample or an array of multiple samples.
"""
pdf(d::AbstractMixtureModel, x::Any)

"""
    logpdf(d::Union{UnivariateMixture, MultivariateMixture}, x)

Evaluate the logarithm of the (mixed) probability density function over `x`.
Here, `x` can be a single sample or an array of multiple samples.
"""
logpdf(d::AbstractMixtureModel, x::Any)

"""
    rand(d::Union{UnivariateMixture, MultivariateDistribution})

Draw a sample from the mixture model `d`.

    rand(d::Union{UnivariateMixture, MultivariateMixture}, n)

Draw `n` samples from `d`.
"""
rand(d::AbstractMixtureModel)

"""
    rand!(d::Union{UnivariateMixture, MultivariateMixture}, r::AbstactArray)

Draw multiple samples from `d` and write them to `r`.
"""
rand!(d::AbstractMixtureModel, r::AbstractArray)

#### Constructors

"""
    MixtureModel(components, [prior])

Construct a mixture model with a vector of `components` and a `prior` probability vector.
If no `prior` is provided then all components will have the same prior probabilities.
"""
MixtureModel(components::Vector{C}) where {C<:Distribution} =
    MixtureModel(components, Categorical(length(components)))

"""
    MixtureModel(C, params, [prior])

Construct a mixture model with component type ``C``, a vector of parameters for constructing
the components given by ``params``, and a prior probability vector.
If no `prior` is provided then all components will have the same prior probabilities.
"""
function MixtureModel(::Type{C}, params::AbstractArray) where C<:Distribution
    components = C[_construct_component(C, a) for a in params]
    MixtureModel(components)
end

function MixtureModel(components::Vector{C}, prior::Categorical) where C<:Distribution
    VF = variate_form(C)
    VS = value_support(C)
    MixtureModel{VF,VS,C}(components, prior)
end

MixtureModel(components::Vector{C}, p::Vector{Float64}) where {C<:Distribution} =
    MixtureModel(components, Categorical(p))

_construct_component(::Type{C}, arg) where {C<:Distribution} = C(arg)
_construct_component(::Type{C}, args::Tuple) where {C<:Distribution} = C(args...)

function MixtureModel(::Type{C}, params::AbstractArray, p::Vector{Float64}) where C<:Distribution
    components = C[_construct_component(C, a) for a in params]
    MixtureModel(components, p)
end

#### Basic properties

"""
    length(d::MultivariateMixture)

The length of each sample (only for `Multivariate`).
"""
length(d::MultivariateMixture) = length(d.components[1])
size(d::MatrixvariateMixture) = size(d.components[1])

ncomponents(d::MixtureModel) = length(d.components)
components(d::MixtureModel) = d.components
component(d::MixtureModel, k::Int) = d.components[k]

probs(d::MixtureModel) = probs(d.prior)
params(d::MixtureModel) = ([params(c) for c in d.components], params(d.prior)[1])
partype(d::MixtureModel) = promote_type(partype(d.prior), map(partype, d.components)...)

minimum(d::MixtureModel) = minimum([minimum(dci) for dci in d.components])
maximum(d::MixtureModel) = maximum([maximum(dci) for dci in d.components])

function mean(d::UnivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    m = 0.0
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            m += mean(c) * pi
        end
    end
    return m
end

function mean(d::MultivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    m = zeros(length(d))
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            BLAS.axpy!(pi, mean(c), m)
        end
    end
    return m
end

"""
    var(d::UnivariateMixture)

Compute the overall variance (only for ``UnivariateMixture``).
"""
function var(d::UnivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    means = Vector{Float64}(undef, K)
    m = 0.0
    v = 0.0
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            ci = component(d, i)
            means[i] = mi = mean(ci)
            m += pi * mi
            v += pi * var(ci)
        end
    end
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            v += pi * abs2(means[i] - m)
        end
    end
    return v
end

function cov(d::MultivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    m = zeros(length(d))
    md = zeros(length(d))
    V = zeros(length(d),length(d))

    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            BLAS.axpy!(pi, mean(c), m)
            BLAS.axpy!(pi, cov(c), V)
        end
    end
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            # todo: use more in-place operations
            md = mean(c) - m
            BLAS.axpy!(pi, md*md', V)
        end
    end
    return V
end

#### show

function show(io::IO, d::MixtureModel)
    K = ncomponents(d)
    pr = probs(d)
    println(io, "MixtureModel{$(component_type(d))}(K = $K)")
    Ks = min(K, 8)
    for i = 1:Ks
        @printf(io, "components[%d] (prior = %.4f): ", i, pr[i])
        println(io, component(d, i))
    end
    if Ks < K
        println(io, "The rest are omitted ...")
    end
end

#### Evaluation

function insupport(d::AbstractMixtureModel, x::AbstractVector)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0 && insupport(component(d, i), x)
            return true
        end
    end
    return false
end

function _cdf(d::UnivariateMixture, x::Real)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    r = 0.0
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            c = component(d, i)
            r += pi * cdf(c, x)
        end
    end
    return r
end

cdf(d::UnivariateMixture{Continuous}, x::Float64) = _cdf(d, x)
cdf(d::UnivariateMixture{Discrete}, x::Int) = _cdf(d, x)

function _mixpdf1(d::AbstractMixtureModel, x)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    v = 0.0
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            c = component(d, i)
            v += pdf(c, x) * pi
        end
    end
    return v
end

function _mixpdf!(r::AbstractArray, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    fill!(r, 0.0)
    t = Array{Float64}(undef, size(r))
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            if d isa UnivariateMixture
                t .= pdf.(component(d, i), x)
            else
                pdf!(t, component(d, i), x)
            end
            BLAS.axpy!(pi, t, r)
        end
    end
    return r
end

function _mixlogpdf1(d::AbstractMixtureModel, x)
    # using the formula below for numerical stability
    #
    # logpdf(d, x) = log(sum_i pri[i] * pdf(cs[i], x))
    #              = log(sum_i pri[i] * exp(logpdf(cs[i], x)))
    #              = log(sum_i exp(logpri[i] + logpdf(cs[i], x)))
    #              = m + log(sum_i exp(logpri[i] + logpdf(cs[i], x) - m))
    #
    #  m is chosen to be the maximum of logpri[i] + logpdf(cs[i], x)
    #  such that the argument of exp is in a reasonable range
    #

    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K

    lp = Vector{Float64}(undef, K)
    m = -Inf   # m <- the maximum of log(p(cs[i], x)) + log(pri[i])
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            # lp[i] <- log(p(cs[i], x)) + log(pri[i])
            lp_i = logpdf(component(d, i), x) + log(pi)
            @inbounds lp[i] = lp_i
            if lp_i > m
                m = lp_i
            end
        end
    end
    v = 0.0
    @inbounds for i = 1:K
        if p[i] > 0.0
            v += exp(lp[i] - m)
        end
    end
    return m + log(v)
end

function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    n = length(r)
    Lp = Matrix{Float64}(undef, n, K)
    m = fill(-Inf, n)
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            lpri = log(pi)
            lp_i = view(Lp, :, i)
            # compute logpdf in batch and store
            if d isa UnivariateMixture
                lp_i .= logpdf.(component(d, i), x)
            else
                logpdf!(lp_i, component(d, i), x)
            end

            # in the mean time, add log(prior) to lp and
            # update the maximum for each sample
            for j = 1:n
                lp_i[j] += lpri
                if lp_i[j] > m[j]
                    m[j] = lp_i[j]
                end
            end
        end
    end

    fill!(r, 0.0)
    @inbounds for i = 1:K
        if p[i] > 0.0
            lp_i = view(Lp, :, i)
            for j = 1:n
                r[j] += exp(lp_i[j] - m[j])
            end
        end
    end

    @inbounds for j = 1:n
        r[j] = log(r[j]) + m[j]
    end
    return r
end

pdf(d::UnivariateMixture{Continuous}, x::Real) = _mixpdf1(d, x)
pdf(d::UnivariateMixture{Discrete}, x::Int) = _mixpdf1(d, x)
logpdf(d::UnivariateMixture{Continuous}, x::Real) = _mixlogpdf1(d, x)
logpdf(d::UnivariateMixture{Discrete}, x::Int) = _mixlogpdf1(d, x)

_pdf!(r::AbstractArray, d::UnivariateMixture{Discrete}, x::UnitRange) = _mixpdf!(r, d, x)
_pdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixpdf!(r, d, x)
_logpdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixlogpdf!(r, d, x)

_pdf(d::MultivariateMixture, x::AbstractVector) = _mixpdf1(d, x)
_logpdf(d::MultivariateMixture, x::AbstractVector) = _mixlogpdf1(d, x)
_pdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixpdf!(r, d, x)
_lodpdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixlogpdf!(r, d, x)

## component-wise pdf and logpdf

function _cwise_pdf1!(r::AbstractVector, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    length(r) == K || error("The length of r should match the number of components.")
    for i = 1:K
        r[i] = pdf(component(d, i), x)
    end
    r
end

function _cwise_logpdf1!(r::AbstractVector, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    length(r) == K || error("The length of r should match the number of components.")
    for i = 1:K
        r[i] = logpdf(component(d, i), x)
    end
    r
end

function _cwise_pdf!(r::AbstractMatrix, d::AbstractMixtureModel, X)
    K = ncomponents(d)
    n = size(X, ndims(X))
    size(r) == (n, K) || error("The size of r is incorrect.")
    for i = 1:K
        if d isa UnivariateMixture
            view(r,:,i) .= pdf.(Ref(component(d, i)), X)
        else
            pdf!(view(r,:,i),component(d, i), X)
        end
    end
    r
end

function _cwise_logpdf!(r::AbstractMatrix, d::AbstractMixtureModel, X)
    K = ncomponents(d)
    n = size(X, ndims(X))
    size(r) == (n, K) || error("The size of r is incorrect.")
    for i = 1:K
        if d isa UnivariateMixture
            view(r,:,i) .= logpdf.(Ref(component(d, i)), X)
        else
            logpdf!(view(r,:,i), component(d, i), X)
        end
    end
    r
end

componentwise_pdf!(r::AbstractVector, d::UnivariateMixture, x::Real) = _cwise_pdf1!(r, d, x)
componentwise_pdf!(r::AbstractVector, d::MultivariateMixture, x::AbstractVector) = _cwise_pdf1!(r, d, x)
componentwise_pdf!(r::AbstractMatrix, d::UnivariateMixture, x::AbstractVector) = _cwise_pdf!(r, d, x)
componentwise_pdf!(r::AbstractMatrix, d::MultivariateMixture, x::AbstractMatrix) = _cwise_pdf!(r, d, x)

componentwise_logpdf!(r::AbstractVector, d::UnivariateMixture, x::Real) = _cwise_logpdf1!(r, d, x)
componentwise_logpdf!(r::AbstractVector, d::MultivariateMixture, x::AbstractVector) = _cwise_logpdf1!(r, d, x)
componentwise_logpdf!(r::AbstractMatrix, d::UnivariateMixture, x::AbstractVector) = _cwise_logpdf!(r, d, x)
componentwise_logpdf!(r::AbstractMatrix, d::MultivariateMixture, x::AbstractMatrix) = _cwise_logpdf!(r, d, x)

componentwise_pdf(d::UnivariateMixture, x::Real) = componentwise_pdf!(Vector{Float64}(undef, ncomponents(d)), d, x)
componentwise_pdf(d::UnivariateMixture, x::AbstractVector) = componentwise_pdf!(Matrix{Float64}(undef, length(x), ncomponents(d)), d, x)
componentwise_pdf(d::MultivariateMixture, x::AbstractVector) = componentwise_pdf!(Vector{Float64}(undef, ncomponents(d)), d, x)
componentwise_pdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_pdf!(Matrix{Float64}(undef, size(x,2), ncomponents(d)), d, x)

componentwise_logpdf(d::UnivariateMixture, x::Real) = componentwise_logpdf!(Vector{Float64}(undef, ncomponents(d)), d, x)
componentwise_logpdf(d::UnivariateMixture, x::AbstractVector) = componentwise_logpdf!(Matrix{Float64}(undef, length(x), ncomponents(d)), d, x)
componentwise_logpdf(d::MultivariateMixture, x::AbstractVector) = componentwise_logpdf!(Vector{Float64}(undef, ncomponents(d)), d, x)
componentwise_logpdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_logpdf!(Matrix{Float64}(undef, size(x,2), ncomponents(d)), d, x)

## Sampling

struct MixtureSampler{VF,VS,Sampler} <: Sampleable{VF,VS}
    csamplers::Vector{Sampler}
    psampler::AliasTable
end

function MixtureSampler(d::MixtureModel{VF,VS}) where {VF,VS}
    csamplers = map(sampler, d.components)
    psampler = sampler(d.prior)
    MixtureSampler{VF,VS,eltype(csamplers)}(csamplers, psampler)
end

rand(rng::AbstractRNG, s::MixtureSampler{Univariate}) =
    rand(rng, s.csamplers[rand(rng, s.psampler)])
rand(rng::AbstractRNG, d::MixtureModel{Univariate}) =
    rand(rng, component(d, rand(rng, d.prior)))

# multivariate mixture sampler for a vector
_rand!(rng::AbstractRNG, s::MixtureSampler{Multivariate}, x::AbstractVector) =
    _rand!(rng, s.csamplers[rand(rng, s.psampler)], x)
_rand!(rng::AbstractRNG, s::MixtureModel{Multivariate}, x::AbstractVector) =
    _rand!(rng, sampler(s), x)

sampler(d::MixtureModel) = MixtureSampler(d)


text:
include("preprocessing.jl")
include("postprocessing.jl")


text:
"""
    AbstractMod{T}

Abstract type for all Modueles of type T

"""
abstract type AbstractMod{T} end

"""
    Mod{T} <: AbstractMod{T}

Structure to store one specific modules of type{T}

"""

# problem input as array in unique gives an error if dims = 2 not added... ne
struct Mod{T} <: AbstractMod{T}
    m::T
end

# some base base function for the new type
"""
    isless(mod1::Mod, mod2::Mod)

Extationtions of ``Base.isless`` for modules of type Mod
for more infromation [see](https://docs.julialang.org/en/v1/base/base/)
"""
Base.isless(mod1::Mod, mod2::Mod) = isless(mod1.m,mod2.m)

"""
    isequal(mod1::Mod, mod2::Mod)

Extationtion of ``Base.isequal`` for modules of type ``Mod``
for more infromation [see](https://docs.julialang.org/en/v1/base/base/)
"""
Base.isequal(mod1::Mod, mod2::Mod) = isequal(mod1.m,mod2.m)

"""
    Base.length(mod1::Mod)

Extationtions of ``Base.length`` for modules of type ``Mod``
"""

Base.length(mod1::Mod) = 1

"""
    Group_Mod{N <: Mod{T} where T} <: AbstractMod{N}

Structure to group multiple modules.
The input values are first filtered to prevent duplicated modules.
Afterwards, the modules are sorted to give consistent results even if modules are load in a different order.
"""

struct Group_Mod{N <: Mod{T} where T} <: AbstractMod{N}
    m::Array{N}
    Group_Mod(m) =  m |> Set |> collect |> sort |> (y -> new{eltype(m)}(y))
end

"""
    group_mod(input::Array{T} where T)

A Function to facilitate the input of multiple modules. It returns a "Group_mod" structure.
The input is an array containing the data that needs to be transformed to an ``Mod`` and grouped afterwards. See
[`Group_mod`](@ref)

```jldoctest
Mods = [:a,:b,:c]
grouped_mods = group_mod(Mods)
isa(grouped_mods,Group_Mod)

# output

true
```
"""
group_mod(input::Array{T} where T) = Group_Mod([Mod(newmod) for newmod in input])

"""
    Base.iterate(Group::Group_Mod,state =1 )

Extationtions of ``Base.iterate`` for modules of type ``Group_Mod``
for more infromation [see](https://docs.julialang.org/en/v1/manual/interfaces/)
"""
function Base.iterate(Group::Group_Mod,state =1 )
    if  state <= length(Group.m)
        mod = Group.m[state]
        state += 1
        return (mod,state)
    else
        return
    end
end

"""
    Base.length(m::Group_Mod)

Return the number of `Mod` in `m`
"""
Base.length(m::Group_Mod) = length(m.m)

"""
    Base.eltype(m::Group_Mod)

Return the type of the iterator
"""
Base.eltype(::Group_Mod{T}) where {T} = T