text:
#==========================================================================================#

# INFLUENCE FUNCTION

function influence(obj::MLE, w::AbstractWeights)
    s = score(obj)
    j = jacobian(obj, w)
    return - s / j'
end

function influence(obj::TwoStageModel, w::AbstractWeights)
    c = influence(first_stage(obj), w) * crossjacobian(obj, w)'
    s = score(obj)
    j = jacobian(obj, w)
    return - (s + c) / j'
end

#==========================================================================================#

# SET COVARIANCE MATRIX

_vcov!(obj::ParModel, V::Matrix) = (switch_stage(obj).V = V)

#==========================================================================================#

# COVARIANCE MATRIX FOR HOMOSCEDASTIC MLE

function _vcov!(obj::MLE, corr::Homoscedastic, w::AbstractWeights)
    if getcorr(obj).method == "OIM"
        obj.V = - inv(jacobian(obj, w))
    elseif getcorr(obj).method == "OPG"
        g     = score(obj)
        obj.V = inv(crossprod(g, w))
    else
        throw("unknown method")
    end
end

#==========================================================================================#

# COVARIANCE MATRIX FOR MLE AND TWO-STAGE MODELS

for m in (MLE, TwoStageModel)
    @eval begin

        function _vcov!(obj::$m, corr::Heteroscedastic, w::UnitWeights)
            ψ = influence(obj, w)
            V = crossprod(ψ)
            adjfactor!(V, obj, corr)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Heteroscedastic, w::FrequencyWeights)
            ψ = influence(obj, w)
            V = crossprod(ψ, w)
            adjfactor!(V, obj, corr)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Heteroscedastic, w::AbstractWeights)
            ψ = influence(obj, w) ; lmul!(Diagonal(w), ψ)
            V = crossprod(ψ)
            adjfactor!(V, obj, corr)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Clustered, w::UnitWeights)
            ψ = corr.mat * influence(obj, w)
            V = crossprod(ψ)
            adjfactor!(V, obj, corr)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::Clustered, w::AbstractWeights)
            ψ = corr.mat * influence(obj, w) ; lmul!(Diagonal(w), ψ)
            V = crossprod(ψ)
            adjfactor!(V, obj, corr)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::CrossCorrelated, w::UnitWeights)
            ψ = influence(obj, w)
            V = crossprod(ψ, corr.mat)
            adjfactor!(V, obj, corr)
            _vcov!(obj, V)
        end

        function _vcov!(obj::$m, corr::CrossCorrelated, w::AbstractWeights)
            ψ = influence(obj, w) ; lmul!(Diagonal(w), ψ)
            V = crossprod(ψ, corr.mat)
            adjfactor!(V, obj, corr)
            _vcov!(obj, V)
        end
    end
end

#==========================================================================================#

# GMM WEIGHT MATRIX

function wmatrix(obj::GMM, corr::Heteroscedastic, ::UnitWeights)
    s = score(obj)
    return crossprod(s)
end

function wmatrix(obj::GMM, corr::Heteroscedastic, w::FrequencyWeights)
    s = score(obj)
    return crossprod(s, w)
end

function wmatrix(obj::GMM, corr::Heteroscedastic, w::AbstractWeights)
    s = score(obj) ; lmul!(Diagonal(w), s)
    return crossprod(s)
end

function wmatrix(obj::GMM, corr::Clustered, ::UnitWeights)
    s = corr.mat * score(obj)
    return crossprod(s)
end

function wmatrix(obj::GMM, corr::Clustered, w::AbstractWeights)
    s = corr.mat * score(obj) ; lmul!(Diagonal(w), s)
    return crossprod(s)
end

function wmatrix(obj::GMM, corr::CrossCorrelated, ::UnitWeights)
    s = score(obj)
    return crossprod(s, corr.mat)
end

function wmatrix(obj::GMM, corr::CrossCorrelated, w::AbstractWeights)
    s = score(obj) ; lmul!(Diagonal(w), s)
    return crossprod(s, corr.mat)
end

# COVARIANCE MATRIX FOR GMM

function _vcov!(obj::GMM, corr::CorrStructure, w::AbstractWeights)

    j = jacobian(obj, w)
    S = wmatrix(obj, corr, w)

    if obj.method == "Method of moments"
        V = j \ (S / j')
    elseif obj.method == "One-step GMM"
        ψ = j / (j' * j)
        V = crossprod(ψ, S)
    elseif obj.method == "Two-step GMM"
        ω = obj.W \ j
        ψ = ω / (j' * ω)
        V = crossprod(ψ, S)
    elseif obj.method == "Optimal GMM"
        ω = S \ j
        V = inv(j' * ω)
    end

    adjfactor!(V, obj, corr)
    _vcov!(obj, V)
end


text:
function FieldLnsCoefficients(arg0::jint, arg1::jint, arg2::Vector{Vector{RealFieldElement}}, arg3::TreeMap, arg4::RealFieldElement, arg5::Field)
    return FieldLnsCoefficients((jint, jint, Vector{Vector{RealFieldElement}}, TreeMap, RealFieldElement, Field), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_lns(obj::FieldLnsCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getLns", RealFieldElement, (jint, jint), arg0, arg1)
end

function getd_lnsd_gamma(obj::FieldLnsCoefficients, arg0::jint, arg1::jint)
    return jcall(obj, "getdLnsdGamma", RealFieldElement, (jint, jint), arg0, arg1)
end

ext: jl


text:
module TestDoctest

using Documenter
using Test
using Transducers

@testset "doctest" begin
    if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
        @info "Skipping doctests on PkgEval."
        return
    elseif VERSION >= v"1.5-"
        @info "Skipping doctests on Julia $VERSION."
        return
    end
    include("../docs/utils.jl")
    transducers_rm_examples()
    doctest(Transducers; manual=true)
end

end  # module


text:
function capitalize(cur_string::AbstractString)
  cur_string = join(map(x -> ucfirst(x), split(cur_string, "/")), "/")

  cur_string
end


text:
module Operators

using StaticArrays: maximum
import ..slab
import ..column
import ..DataLayouts: DataLayouts, Data2D, DataSlab2D
import ..Geometry:
    Geometry, Cartesian12Vector, Covariant12Vector, Contravariant12Vector, ⊗
import ..Spaces: Spaces, Quadratures, AbstractSpace
import ..Topologies
import ..Fields: Fields, Field
using ..RecursiveApply

import LinearAlgebra
using StaticArrays

include("spectralelement.jl")
include("numericalflux.jl")
include("finitedifference.jl")

end # module


text:
@safetestset "blocking" begin
    using DifferentialEvolutionMCMC, Test, Random, Distributions
    using MCMCChains
    Random.seed!(58122)

    data = rand(Normal(0.0, 1.0), 100)

    function prior_loglike(μ, σ)
        LL = 0.0
        LL += logpdf(Normal(0, 10), μ)
        LL += logpdf(truncated(Cauchy(0, 1), 0, Inf), σ)
        return LL
    end

    function sample_prior()
        μ = rand(Normal(0, 10))
        σ = rand(truncated(Cauchy(0, 1), 0, Inf))
        return [μ,σ]
    end

    bounds = ((-Inf,Inf),(0.0,Inf))

    data = rand(Normal(0.0, 1.0), 1000)

    function loglike(data, μ, σ)
        return sum(logpdf.(Normal(μ, σ), data))
    end

    bounds = ((-Inf,Inf),(0.0,Inf))

    names = (:μ,:σ)

    blocks = [[true,false],[false,true]]

    blocks = as_union(blocks)

    blocking_on = x -> true

    model = DEModel(;
        sample_prior,
        prior_loglike,
        loglike,
        data,
        names,
    )

    de = DE(;
        sample_prior,
        bounds,
        burnin = 1000,
        Np = 6,
        blocking_on,
        blocks
    )

    n_iter = 2000
    chains = sample(model, de, n_iter, progress=true)

    means = describe(chains)[1][:,:mean]
    rhat = describe(chains)[1][:,:rhat]

    @test means[1] ≈ 0.0 atol = .1
    @test means[2] ≈ 1.0 atol = .1
    @test rhat[1] ≈ 1.0 atol = .01
    @test rhat[2] ≈ 1.0 atol = .01

    chains = sample(model, de, MCMCThreads(), n_iter, progress=true)

    means = describe(chains)[1][:,:mean]
    rhat = describe(chains)[1][:,:rhat]

    @test means[1] ≈ 0.0 atol = .1
    @test means[2] ≈ 1.0 atol = .1
    @test rhat[1] ≈ 1.0 atol = .01
    @test rhat[2] ≈ 1.0 atol = .01
end


text:
using NeuralDataScience
using Documenter

DocMeta.setdocmeta!(NeuralDataScience, :DocTestSetup, :(using NeuralDataScience); recursive=true)

makedocs(;
    modules=[NeuralDataScience],
    authors="Martin Wimpff",
    repo="https://github.com/martinwimpff/NeuralDataScience.jl/blob/{commit}{path}#{line}",
    sitename="NeuralDataScience.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://martinwimpff.github.io/NeuralDataScience.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/martinwimpff/NeuralDataScience.jl",
    devbranch="main",
)


text:
##############################################################################
# This code generates all the figures in the portfolio slides.
#
# This material was developed as part of the course EE103 taught at Stanford
# University by Professor Stephen Boyd. For details, see the course
# website at <web.stanford.edu/class/ee103/>.
##############################################################################

using PyPlot;
using LinearLeastSquares;

##############################################################################
# First load the stock data and form the returns matrix
##############################################################################
include("portfolio_optimization_data.jl");
T, n = size(p)# Number of periods and number of assets

##############################################################################
#
# Single Asset Investment
#
##############################################################################
assets = ["KO", "BP"]

#Plot prices of the stocks over the last 10 years
ind1 = find(stockNames.==assets[1]);
ind2 = find(stockNames.==assets[2]);
figure();
ax = gca();
ax[:set_ylim]([0,70]);
ax[:set_xlim]([0,T]);
plot(p[:,ind1], label = assets[1])
plot(p[:,ind2], label = assets[2]);
xlabel("Days");
ylabel("Prices");
legend();
savefig("stock_prices.eps")

#Plot price changes over a few weeks
figure();
ax = gca();
ax[:set_ylim]([0,70]);
ax[:set_xlim]([1600,1650]);
plot(p[:,ind1], label = assets[1])
plot(p[:,ind2], label = assets[2]);
xlabel("Days");
ylabel("Prices");
legend();
savefig("stock_prices_zoomed.eps")

#Plot return of assets over the same period (a few weeks)
figure();
ax = gca();
ax[:set_ylim]([-0.2,0.2]);
ax[:set_xlim]([1600,1650]);
plot(R[:,ind1], label = assets[1])
plot(R[:,ind2], label = assets[2]);
xlabel("Days");
ylabel("Returns");
legend();
savefig("stock_return_zoomed.eps")

#Generate cumulative returns for two assets
#with an initial budget
assets = ["KO", "MSFT"];
budget = 1e4;
ind1 = find(stockNames.==assets[1]);
ind2 = find(stockNames.==assets[2]);
h1 = (budget/p[1, ind1][1])*p[:, ind1];
h2 = (budget/p[1, ind2][1])*p[:, ind2];

#Plot the cumulative returns
figure();
ax = gca();
ax[:set_ylim]([5000,25000]);
ax[:set_xlim]([0,T]);
plot(1:T,h1, label = assets[1]);
plot(1:T, h2, label = assets[2]);
xlabel("Days");
ylabel("Value");
legend();
savefig("cumulative_return.eps")

#Plot risk returns for all assets
figure();
ax = gca();
ax[:set_ylim]([0,25]);
ax[:set_xlim]([0,60]);
assets = stockNames; #risk return string for all assets
for i = 1:length(assets)
    k = find(stockNames.==assets[i]);
    r = R[:, k];
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    scatter(annualRisk, annualReturn);
end

#Annotate some of the stocks on the scatter plot
annotateStocks = ["USDOLLAR", "SBUX", "GS", "MMM", "BRCM"];
for(i,txt) in enumerate(annotateStocks)
	ind = find(stockNames .== txt);
	ax[:annotate](txt, (risks[ind], returns[ind]));
end
xlabel("Annualized Risk");
ylabel("Annualized Return");
savefig("asset_risk_return.eps");

##############################################################################
# Define a function to perform
# constant weight portfolio optimization
# Returns vector w minimizing norm(R*w - rho*ones)^2
# subject to sum(w) = 1, mean(R*w) = rho
##############################################################################
function optimalPortfolio(R, rho)
	w = Variable(size(R,2));
	constraints = [sum(w) == 1, mean(R*w) == rho];
	minimize!(sum_squares(R*w - rho), constraints);
	return evaluate(w);
end
# we use this function for the rest
# of the figure generation

##############################################################################
#
# Portfolio Investment
#
##############################################################################

#risk return plot of various assets

#Generate a "good" constant value portfolio
rho = 0.18/250;
w_good= convert(Array,full(optimalPortfolio(R, rho)));

#Generate a bad portfolio
#first, get indicies corresponding to some bad assets
assets = ["HPQ", "KO", "MSFT", "BAC"];
ks = Int64[];
for i = 1:length(assets)
    push!(ks, find(stockNames .== assets[i])[1] );
end
#invest in "bad" stocks and short in "good" stocks
w_bad = zeros(n, 1);
w_bad[ks] = [8000; -4000; 10000; -4000];

#generate a uniform weight portfolio
w_uniform = (ones(1,n)./n)';

#Create the holdings vector
#for each portfolio,
#uniform, good, bad
BUDGET = 10000;
H = BUDGET.*[w_uniform w_good w_bad];
labels = ["uniform", "good", "bad"];

figure();
ax = gca();
ax[:set_ylim]([0,25]);
ax[:set_xlim]([0,60]);

#generate the risk/return points of these portfolios
# alongside all single value investment portfolios
assets = stockNames; #risk return string for all assets
for i = 1:length(assets)
    k = find(stockNames.==assets[i]);
    r = R[:, k];
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    scatter(annualRisk, annualReturn)
end

#generate risk/return points for uniform/good/bad
for i = 1:size(H, 2)
	txt = labels[i];
    h = H[:, i];
    annualRisk = sqrt(250)*100*std(R*h)/sum(h);
    annualReturn = 250*100*mean(R*h)/sum(h);
    ax[:annotate](txt, (annualRisk, annualReturn));
    scatter(annualRisk, annualReturn, c = "red")
end

xlabel("Annualized Risk");
ylabel("Annualized Return");
savefig("portfolio_risk_return.eps")

#### Constant weight portfolio with reinvestment
assets = ["KO", "MSFT"];
ks = Int64[];
for i = 1:length(assets)
	k = find(stockNames.==assets[i])[1];
	push!(ks,k);
end

w = zeros(n, 1);
w[ks] = [0.5; 0.5]; # uniform across KO and MSFT
V = Float64[];
push!(V,10000);
for t = 1:T-1
    h = V[end]*w;
    push!(V, V[end] + (R[t, :]*h)[1]);
end

figure();
plot(V, label = "uniform portfolio");
plot(p[:, ks[1]]*10000/p[1, ks[1]], label = assets[1], ls="--");
plot(p[:, ks[2]]*10000/p[1, ks[2]], label = assets[2],ls="--");

ax = gca();
ax[:set_ylim]([0,3e4]);
ax[:set_xlim]([0,T]);
xlabel("Days");
ylabel("Value");
legend();
savefig("const_weight_cumulative_return.eps")

#Heavily leveraged portfolio with reinvestment (going bust)
assets = ["KO", "MSFT"];
ks = Int64[];
for i = 1:length(assets)
	k = find(stockNames.==assets[i])[1];
	push!(ks,k);
end

w = zeros(n, 1);
w[ks] = [-3; 4];
V = Float64[];
push!(V,10000);
for t = 1:T-1
    h = V[end]*w;
    push!(V, V[end] + (R[t, :]*h)[1]);
end

#after this point we stop investing
bust_index = find(V .< 1000)[1];

figure();
plot(V[1:bust_index]);
plot(p[:, ks[1]]*10000/p[1, ks[1]], label = assets[1], ls="--");
plot(p[:, ks[2]]*10000/p[1, ks[2]], label = assets[2],ls="--");
plot(bust_index:T, 1000*ones(T-bust_index+1), label = "leveraged portfolio", c = "red");

ax = gca();
ax[:set_ylim]([0,3e4]);
ax[:set_xlim]([0,T]);
xlabel("Days");
ylabel("Value");
legend();
savefig("bust_cumulative_return.eps")

# Reinvestment or not
h = 10000;
V1 = Float64[];  #without reinvestment
V2 = Float64[]; #with reinvestment
push!(V1, h);
push!(V2, h);

k = find(stockNames .=="COST")[1];

for t = 1:T-1
	push!(V1, V1[t] + R[t,k]*h);
	push!(V2, V2[t]*(1 + R[t,k]));
end
figure();
plot(1:T, V1, label = "no reinvestment");
plot(1:T, V2, label = "reinvestment");

ax = gca();
ax[:set_ylim]([0,5e4]);
ax[:set_xlim]([0,T]);
xlabel("Days");
ylabel("Value");
legend();
savefig("reinvestment_comparison.eps")

##############################################################################
#
# Portfolio Optimization
#
##############################################################################
## cumulative value plot for some optimal portfolios
w20 = optimalPortfolio(R, 0.20/250);
w25 = optimalPortfolio(R, 0.25/250);
V20, V25 = Float64[], Float64[];
push!(V20, 10000);
push!(V25, 10000);
for t = 1:(T-1)
   push!(V20, (V20[t]*(1 + R[t,:]*w20))[1]    );
   push!(V25, (V25[t]*(1 + R[t,:]*w25))[1]   );
end

figure();
plot(1:T, V20, label = "optimal portfolio, rho = 0.20/250");
plot(1:T, V25, c = "g", label = "optimal portfolio, rho=0.25/250");
plot(1:T, p[:, 3]*10000/p[1, 3], c ="r", label = "individual assets");
plot(1:T, p[:, 5]*10000/p[1, 5], c = "r");
plot(1:T, p[:, 7]*10000/p[1, 7], c = "r");
plot(1:T, p[:, 8]*10000/p[1, 8], c = "r");

ax = gca();
ax[:set_ylim]([5e3,5e5]);
ax[:set_xlim]([0,T]);
ax[:set_yscale]("log");
xlabel("Days");
ylabel("Value");
legend();
savefig("optimal_cumulative_value.eps")

##### Risk-return plot of assets
min_return = mean(R[:, end]); # return of the risk-free asset
rhos = linspace(min_return, 0.25/250, 11);
risks, returns = Float64[], Float64[];
figure();
for rho in rhos
    w = optimalPortfolio(R, rho);
    r = R*w;
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    push!(risks, annualRisk);
    push!(returns, annualReturn);
    scatter(annualRisk, annualReturn, s= 100, c = "r");
end

plot(risks, returns, c = "r");
for i = 1:length(stockNames)
    r = R[:, i];
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    scatter(annualRisk, annualReturn, s = 100, c = "b")
end

ax = gca();
xlabel("Annualized Risk");
ylabel("Annualized Return");
ax[:set_ylim]([0,25.1]);
ax[:set_xlim]([0,60]);
savefig("optimal_risk_return_curve.eps")

## Risk-return plot of assets (train-test), when BA works reasonably
min_return = mean(R[:, end]); # return of the risk-free asset
rhos = linspace(min_return*2, 0.25/250, 7);

split_day = 2000;
T_train = 900;
T_test = 200;
R_train = R[split_day-T_train:split_day-1, :];
R_test = R[split_day:split_day+T_test-1, :];

trainResult, testResult = zeros(length(rhos), 2), zeros(length(rhos),2);
for (i, rho) in enumerate(rhos)
    w = optimalPortfolio(R_train, rho);
    r = R_train*w;
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    trainResult[i,:] = [annualRisk annualReturn];
    r = R_test*w;
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    testResult[i,:] = [annualRisk annualReturn];
end

figure();
scatter(trainResult[:, 1], trainResult[:, 2], 100, c = "r", label = "Train");
scatter(testResult[:, 1], testResult[:, 2], 100, c = "b", label = "Test");
xlabel("Annualized Risk");
ylabel("Annualized Return");
ax = gca();
ax[:set_ylim]([0,27]);
ax[:set_xlim]([0,16]);
legend();
savefig("train_test_comparison.eps")

# Show the interval
figure();
p_norm = p;
for i in 1:n
    p_norm[:, i] = p_norm[:, i]*BUDGET/p_norm[1, i];
end
plot(p_norm[:, 1:3:n], c= "r");
axvline(x = split_day)
axvline(x = split_day - T_train);
axvline(x = split_day + T_test);
ax = gca();
ax[:annotate]("Train", [(split_day-T_train/2);4e4])
ax[:annotate]("Test", [(split_day + T_test/4);4e4])
ax[:set_ylim]([0,5e4]);
ax[:set_xlim]([0,T]);
legend();
savefig("train_test_comparison_interval.eps")

## Risk-return plot of assets (train-test), when BA doesn't work out
min_return = mean(R[:, end]); # return of the risk-free asset
rhos = linspace(min_return*2, 0.20/250, 7);

split_day = 1400;
T_train = 900;
T_test = 200;
R_train = R[split_day-T_train:split_day-1, :];
R_test = R[split_day:split_day+T_test-1, :];

trainResult, testResult = zeros(length(rhos), 2), zeros(length(rhos),2);
for (i,rho) in enumerate(rhos)
    # minimizing norm(R_train*w - rho*ones)
    w = optimalPortfolio(R_train, rho);
    r = R_train*w;
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    trainResult[i,:] = [annualRisk annualReturn];
    r = R_test*w;
    annualRisk = sqrt(250)*100*std(r);
    annualReturn = 250*100*mean(r);
    testResult[i,:] = [annualRisk annualReturn];
end

figure();
scatter(trainResult[:, 1], trainResult[:, 2], 100, c = "r", label = "Train");
scatter(testResult[:, 1], testResult[:, 2], 100, c = "b", label = "Test");
xlabel("Annualized Risk");
ylabel("Annualized Return");
ax = gca();
ax[:set_ylim]([-20,27]);
ax[:set_xlim]([0,16]);
legend();
savefig("train_test_comparison_bad.eps")

## Show the interval
figure();
p_norm = p;
for i in 1:n
    p_norm[:, i] = p_norm[:, i]*BUDGET/p_norm[1, i];
end
plot(p_norm[:, 1:3:n], c= "r");
axvline(x = split_day)
axvline(x = split_day - T_train);
axvline(x = split_day + T_test);
ax = gca();
ax[:annotate]("Train", [(split_day-T_train/2);4e4])
ax[:annotate]("Test", [(split_day + T_test/4);4e4])
ax[:set_ylim]([0,5e4]);
ax[:set_xlim]([0,T]);
legend();
savefig("train_test_comparison_bad_interval.eps")

##############################################################################
#
# Daily rolling portfolio
#
##############################################################################
bad_periods = 1200;
p = p[bad_periods:end, :]; # remove the bad periods
(T, n) = size(p);
R = diff(p) ./ p[1:end-1, :];

#rolling portfolio simulation
rhos = (1/(250*100)) * [5, 10, 15];
# Number of periods to use
L = 400;
Vs = zeros(T-L-1,length(rhos));
for (j,rho) in enumerate(rhos)
	h = 10000; # Initial budget
	V = zeros(T-1-L);
	V[1] = h;
	for (i,t) in enumerate(L+1:T-1)
		println(i)
		w = optimalPortfolio(R[t-L:t-1,:], rho);
		r = (R[t, :]*w)[1];
		h = (1+r)*h;
		V[i] = h;
	end
	println(size(Vs));
	println(size(V));
	Vs[:, j] = V;
end

# Plotting
figure();
startInd = bad_periods+L+1;
endInd = bad_periods + T-1;
plot(startInd:endInd, Vs[:,1],label="rho=0.05/250");
plot(startInd:endInd, Vs[:,2],label="rho=0.01/250");
plot(startInd:endInd, Vs[:,3],label="rho=0.15/250");
xlabel("Days");
ylabel("Value");
ax = gca();
ax[:set_ylim]([0.95e4,1.3e4]);
ax[:set_xlim]([bad_periods+L,bad_periods+T]);
legend();
save_figure("rolling_portfolio_daily.eps");

## same rolling portfolio simulation with different update period
rhos = (1/(250*100)) * [5, 10, 15];
# Number of periods to use
L = 400;
K=60; #update every K periods
Vs = zeros(T-L-1,length(rhos));
for (j,rho) in enumerate(rhos)
	h = 10000; # Initial budget
	V = zeros(T-1-L);
	V[1] = h;
	for (i,t) in enumerate(L+1:T-1)
		if(rem(i-1,K) == 0)
			w = optimalPortfolio(R[t-L:t-1,:], rho);
		end
		r = (R[t, :]*w)[1];
		h = (1+r)*h;
		V[i] = h;
	end
	Vs[:, j] = V;
end

# Plotting
figure();
startInd = bad_periods+L+1;
endInd = bad_periods + T-1;
plot(startInd:endInd, Vs[:,1],label="rho=0.05/250");
plot(startInd:endInd, Vs[:,2],label="rho=0.01/250");
plot(startInd:endInd, Vs[:,3],label="rho=0.15/250");
xlabel("Days");
ylabel("Value");
ax = gca();
ax[:set_ylim]([0.95e4,1.3e4]);
ax[:set_xlim]([bad_periods+L,bad_periods+T]);
legend();
save_figure("rolling_portfolio_quarterly.eps");


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

@doc raw"""
	thermal(EP::Model, inputs::Dict, UCommit::Int, Reserves::Int, CapacityReserveMargin::Int)
The thermal module creates decision variables, expressions, and constraints related to thermal power plants e.g. coal, oil or natural gas steam plants, natural gas combined cycle and combustion turbine plants, nuclear, hydrogen combustion etc.
This module uses the following 'helper' functions in separate files: ```thermal_commit()``` for thermal resources subject to unit commitment decisions and constraints (if any) and ```thermal_no_commit()``` for thermal resources not subject to unit commitment (if any).
"""
function thermal(EP::Model, inputs::Dict, UCommit::Int, Reserves::Int, CapacityReserveMargin::Int)

	println("Thermal Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	THERM_COMMIT = inputs["THERM_COMMIT"]
	THERM_NO_COMMIT = inputs["THERM_NO_COMMIT"]
	T = inputs["T"]
	THERM_ALL = inputs["THERM_ALL"]
	dfGen = inputs["dfGen"]

	if !isempty(THERM_COMMIT)
		EP = thermal_commit(EP::Model, inputs::Dict, Reserves::Int)
	end

	if !isempty(THERM_NO_COMMIT)
		EP = thermal_no_commit(EP::Model, inputs::Dict, Reserves::Int)
	end

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

	## CO2 Policy Module Thermal Generation by zone
	@expression(EP, eGenerationByThermAll[z=1:Z, t=1:T], # the unit is GW
		sum(EP[:vP][y,t] for y in intersect(inputs["THERM_ALL"], dfGen[dfGen[!,:Zone].==z,:R_ID]))
	)
	EP[:eGenerationByZone] += eGenerationByThermAll

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
import MPI
using CUDA

precompile(Tuple{typeof(Diffusion.julia_diffusion), UInt64, UInt64, UInt64, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, CuPtr{Cdouble}, MPI.MPI_Comm})

max_stars_repo_name: omlins/libdiffusion

max_issues_repo_name: omlins/libdiffusion

max_forks_repo_name: omlins/libdiffusion


text:
# KnetArray implements GPU arrays.  Important differences from the
# alternative CudaArray: (1) The memory manager minimizes the slow
# cudaMalloc by reusing already allocated but garbage collected
# arrays.  (2) getindex handles ranges as views (with shared memory)
# instead of copies.

# KnetPtr type holds a gpu (dev>=0) or cpu (dev=-1) allocated pointer.

type KnetPtr
    ptr::Ptr{Void}
    len::Int
    dev::Int
    parent
end

# We use the KnetPtrs type to keep track of allocated pointers: We
# need one per size per device.  KnetFree[i+2] will hold a dictionary
# from sizes to KnetPtrs for device i.

type KnetPtrs
    used::Int
    free::Array{Ptr{Void},1}
    KnetPtrs()=new(0,Array(Ptr{Void},0))
end

# When Julia gc reclaims a KnetPtr object, the finalizer does not
# actually release the memory, but inserts it in the KnetFree dict
# keyed by length in bytes so it can be reused.

if !isdefined(:KnetFree)
    KnetFree = [ Dict{Int,KnetPtrs}() for i=1:gpucount()+1 ]
end

function freeKnetPtr(p::KnetPtr)
    ptrs = KnetFree[p.dev+2][p.len]
    push!(ptrs.free,p.ptr)
end

# The KnetPtr constructor tries to avoid actual allocation which is
# slow.  First it tries to find a previously allocated and garbage
# collected pointer in KnetFree[dev+2].  If not available it allocates
# a new one if we have more than 10^8 bytes of free memory.  Otherwise
# it tries running gc() and see if we get a pointer back.  Finally if
# all else fails, it calls knetgc which cleans up all allocated
# KnetPtrs on the current device and starts over.

function KnetPtr(nbytes::Integer)
    dev = gpu()
    ptrs = get!(KnetPtrs,KnetFree[dev+2],nbytes)
    if !isempty(ptrs.free)
        kp = KnetPtr(pop!(ptrs.free),nbytes,dev,nothing)
    elseif dev == -1 || gpufree() > 10^8
        ptr = knetMalloc(nbytes)
        kp = KnetPtr(ptr,nbytes,dev,nothing)
        ptrs.used += 1
    elseif (print("."); gc(); !isempty(ptrs.free))
        kp = KnetPtr(pop!(ptrs.free),nbytes,dev,nothing)
    else
        print((:knetgc,nbytes)); gpuinfo()
        knetgc()
        ptr = knetMalloc(nbytes)
        kp = KnetPtr(ptr,nbytes,dev,nothing)
        ptrs.used += 1
    end
    finalizer(kp, freeKnetPtr)
    return kp
end

# This is used to create a shared pointer.  We need to have the parent field to prevent premature gc.
function KnetPtr(parent::KnetPtr, offs::Integer, len::Integer)
    if len < 0 || offs < 1 || offs+len-1 > parent.len
        throw(BoundsError())
    end
    KnetPtr(parent.ptr+offs-1, len, parent.dev, parent)
end

# If you really want to clean up memory you need to call knetgc:
# Note that this only cleans the current device.

function knetgc()
    dev = gpu()
    gc_enable(false)
    for v in values(KnetFree[dev+2])
        if dev >= 0
            for p in v.free
                @cuda(cudart,cudaFree,(Ptr{Void},),p)
            end
        end
        v.used -= length(v.free)
        empty!(v.free)
    end
    gc_enable(true)
end

function knetMalloc(nbytes::Int)
    if gpu() >= 0
        ptr = Ptr{Void}[C_NULL]
        @cuda(cudart,cudaMalloc,(Ptr{Ptr{Void}},Csize_t),ptr,nbytes)
        return ptr[1]
    else
        convert(Ptr{Void}, pointer(Array(UInt8,nbytes)))
    end
end

meminfo()=[(k,v.used,length(v.free)) for (k,v) in KnetFree[gpu()+2]]

### KnetArray ###

"""

KnetArray is a container for GPU arrays that supports most of the
AbstractArray interface.  Important differences from the alternative
CudaArray are: (1) a custom memory manager that minimizes the number
of calls to the slow cudaMalloc by reusing already allocated but
garbage collected GPU pointers.  (2) a custom getindex that handles
ranges such as `a[5:10]` as views (with memory shared with the
original array) instead of copies.  KnetArrays can be created by
specifying the element type and dimensions or by conversion from
regular Arrays and they can be converted back to regular Arrays (which
involve copying to and from the GPU memory):

    a = KnetArray(Float32,2,3)
    b = KnetArray(zeros(2,3))
    c = Array(b)

"""
type KnetArray{T,N}
    ptr::KnetPtr
    dims::NTuple{N,Int}
end

# Note: I removed <: AbstractArray{T,N} after I painfully discovered
# some inefficient AbstractArray methods inherited unintentionally.
# It is better to define a few extra methods to keep a tighter control
# on what methods exactly get called for KnetArrays.

# Aliases:
typealias KnetMatrix{T} KnetArray{T,2}
typealias KnetVector{T} KnetArray{T,1}
typealias KnetVecOrMat{T} Union{KnetVector{T}, KnetMatrix{T}}

# Constructors:
KnetArray{T,N}(::Type{T}, dims::NTuple{N,Int})=KnetArray{T,N}(KnetPtr(sizeof(T)*prod(dims)), dims)
KnetArray(T::Type, dims::Int...)=KnetArray(T,dims)
KnetArray(T::Type, d::Integer...)=KnetArray(T,convert(Tuple{Vararg{Int}}, d))

# Conversions:
import Base: convert, reshape, unsafe_convert, pointer
# KnetArray <- KnetArray
convert{T,N}(::Type{KnetArray}, x::KnetArray{T,N}) = x
convert{T,N}(::Type{KnetArray{T}}, x::KnetArray{T,N}) = x
convert{T,N}(::Type{KnetArray{T,N}}, x::KnetArray{T,N}) = x
convert{T,N,S}(::Type{KnetArray{T}}, x::KnetArray{S,N}) = convert(KnetArray{T,N}, x)
convert{T,N,S}(::Type{KnetArray{T,N}}, x::KnetArray{S,N}) = convert(KnetArray{T,N},unsafe_copy!(Array(S, size(x)), 1, x, 1, length(x)))
reshape{T}(a::KnetArray{T},dims::Dims)=(if dims==size(a); a; elseif prod(dims)!=length(a); throw(DimensionMismatch()); else; KnetArray{T,length(dims)}(a.ptr,dims); end)
reshape(a::KnetArray, dims::Int...) = reshape(a, dims)
# KnetArray <- AbstractArray
convert{T,N}(::Type{KnetArray}, x::AbstractArray{T,N}) = convert(KnetArray{T,N}, x)
convert{T,N,S}(::Type{KnetArray{T}}, x::AbstractArray{S,N}) = convert(KnetArray{T,N}, x)
convert{T,N,S}(::Type{KnetArray{T,N}}, x::AbstractArray{S,N}) = unsafe_copy!(KnetArray(T, size(x)), 1, convert(Array{T,N},x), 1, length(x))
# Array <- KnetArray
convert{T,N}(::Type{Array}, x::KnetArray{T,N}) = convert(Array{T,N}, x)
convert{T,N,S}(::Type{Array{T}}, x::KnetArray{S,N}) = convert(Array{T,N}, x)
convert{T,N,S}(::Type{Array{T,N}}, x::KnetArray{S,N}) = convert(Array{T,N},unsafe_copy!(Array(S, size(x)), 1, x, 1, length(x)))
# Ptr <- KnetArray
unsafe_convert{T}(::Type{Ptr{T}}, a::KnetArray) = unsafe_convert(Ptr{T}, pointer(a))
pointer{T}(a::KnetArray{T})=convert(Ptr{T}, a.ptr.ptr)
pointer{T}(a::KnetArray{T},i)=convert(Ptr{T}, a.ptr.ptr + (i-1)*sizeof(T))

# AbstractArray interface
import Base: eachindex, elsize, eltype, endof, fill!, first, length, linearindexing, ndims, ones, similar, size, stride, zeros
eachindex(a::KnetArray) = (1:length(a))
elsize{T}(::KnetArray{T}) = sizeof(T)
eltype{T}(::KnetArray{T})=T
endof(a::KnetArray) = length(a)
fill!{T}(a::KnetArray{T},x)=(knetfill!(a,T(x),1,length(a));a)
first(a::KnetArray) = a[1]
length(a::KnetArray)=prod(size(a))
linearindexing(::KnetArray)=LinearFast()
ndims{T,N}(a::KnetArray{T,N})=N
ones{T}(a::KnetArray{T})=fill!(similar(a),one(T))
similar(a::KnetArray, T, dims::Dims)      = KnetArray(T, dims)
similar(a::KnetArray, T, dims::Int...)    = similar(a, T, dims)
similar(a::KnetArray, T)                  = similar(a, T, size(a))
similar{T}(a::KnetArray{T})               = similar(a, T, size(a))
similar{T}(a::KnetArray{T}, dims::Dims)   = similar(a, T, dims)
similar{T}(a::KnetArray{T}, dims::Int...) = similar(a, T, dims)
size(a::KnetArray)=a.dims
size{T,N}(a::KnetArray{T,N},i::Integer)=(if i>N; 1; else; size(a)[i]; end)
stride{T,N}(a::KnetArray{T,N},i::Integer)=(if i>N; length(a); else; s=1; for n=1:(i-1); s*=size(a,n); end; s; end)
strides{T,N}(a::KnetArray{T,N})=ntuple(n->stride(a,n), N)
zeros{T}(a::KnetArray{T})=fill!(similar(a),zero(T))

# Indexing:
import Base: getindex, setindex!

# We will implement indexing ranges as views not copies, if possible (when contiguous).
# For contiguous memory without stride all but the last >1 dimension must be full

# The original getindex(a,i:j...) for AbstractArray copies:
# function _getindex(l::LinearIndexing, A::AbstractArray, I::Union{Real, AbstractArray, Colon}...)
# in abstractarray.jl:487,multidimensional.jl:184.

# which getindex ops does array implement?
# getindex(A::Array, i1::Real)
# getindex(A::Array, i1::Real, i2::Real, I::Real...)
# getindex(A::Array, I::UnitRange{Int})
# getindex(A::Array, c::Colon)
# getindex{T<:Real}(A::Array, I::Range{T})

# Julia #14770
# If I is shorter than ndims(A) but longer than 1 the remaining indices assumed =1
# Also extra 1's at the end of I are ignored

# These two are not sufficient in spite of what the documentation says:
# display goes into an infinite loop!
# getindex{T}(A::KnetArray{T}, i::Int)=unsafe_copy!(T[0], 1, A, i, 1)[1]
# setindex!{T}(A::KnetArray{T}, v, i::Int)=unsafe_copy!(A, i, T[v], 1, 1)

# First deal with the easy cases: integer indices, a Colon or a UnitRange.

function getindex{T}(A::KnetArray{T}, I::Real)
    J = Int(I)
    1 <= J <= length(A) || throw(BoundsError(A,J))
    unsafe_copy!(T[0], 1, A, J, 1)[1]
end

function setindex!{T}(A::KnetArray{T}, v, I::Real)
    J = Int(I)
    1 <= J <= length(A) || throw(BoundsError(A,J))
    unsafe_copy!(A, J, T[v], 1, 1)
end

function getindex{T}(A::KnetArray{T}, I::Real...)
    J = Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        1 <= J[j] <= size(A,j) || throw(BoundsError(A,J))
    end
    i = sub2ind(size(A), J...)
    unsafe_copy!(T[0], 1, A, i, 1)[1]
end

function setindex!{T}(A::KnetArray{T}, v, I::Real...)
    J = Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        1 <= J[j] <= size(A,j) || throw(BoundsError(A,J))
    end
    i = sub2ind(size(A), J...)
    unsafe_copy!(A, i, T[v], 1, 1)
end

function getindex{T}(A::KnetArray{T}, I::UnitRange)
    1 <= first(I) <= last(I) <= length(A) || throw(BoundsError(A,I))
    off = 1+(first(I)-1)*sizeof(T)
    len = length(I)*sizeof(T)
    ptr = KnetPtr(A.ptr, off, len)
    KnetArray{T,1}(ptr, (length(I),))
end

function setindex!{T}(A::KnetArray{T}, v, I::UnitRange)
    1 <= first(I) <= last(I) <= length(A) || throw(BoundsError(A,I))
    if isa(v,Number)
        knetfill!(A,T(v),first(I),length(I))
    elseif (isa(v,KnetArray) || isa(v,Array))
        length(v)==length(I) || throw(DimensionMismatch())
        eltype(v)==T || (v = convert(Array{T},v))
        unsafe_copy!(A,first(I),v,1,length(I))
    else
        throw(MethodError(setindex!, A, v, I))
    end
end

function getindex(A::KnetArray, I::Colon)
    reshape(A,length(A))
end

function setindex!{T}(A::KnetArray{T}, v, I::Colon)
    if isa(v,Number)
        knetfill!(A,T(v),1,length(A))
    elseif (isa(v,KnetArray) || isa(v,Array))
        length(v)==length(A) || throw(DimensionMismatch())
        eltype(v)==T || (v = convert(Array{T},v))
        unsafe_copy!(A,1,v,1,length(A))
    else
        throw(MethodError(setindex!, A, v, I))
    end
end

# TODO: the following getindex, setindex! work for 1 and 2 dimensions only, write general versions.

function getindex{T,N}(A::KnetArray{T,N}, I::Union{Real, UnitRange, Colon}...)
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I...)
    B1 = isa(I[1],Colon) ? size(A,1) : length(I[1])
    B2 = isa(I[2],Colon) ? size(A,2) : length(I[2])
    if ncols == 1
        off = 1+(firstindex-1)*sizeof(T)
        len = nrows*sizeof(T)
        ptr = KnetPtr(A.ptr, off, len)
        KnetArray{T,2}(ptr, (B1,B2))
    else
        B = similar(A, (B1,B2))
        nrows *= sizeof(T); astep *= sizeof(T)
        ccall((:xcopy,libknet8),Void,(Cint,Cint,Ptr{Void},Cint,Ptr{Void},Cint),
              nrows, ncols, pointer(A,firstindex), astep, B, nrows)
        return B
    end
end

function setindex!{T,N}(A::KnetArray{T,N}, B, I::Union{Real, UnitRange, Colon}...)
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I...)
    aptr0 = pointer(A, firstindex)
    if isa(B,Number)
        B = T(B)
        if ncols == 1
            if T <: Float32;    ccall((:fill_32,libknet8),Void,(Cint,Cfloat, Ptr{Cfloat}), nelts,B,aptr0)
            elseif T<: Float64; ccall((:fill_64,libknet8),Void,(Cint,Cdouble,Ptr{Cdouble}),nelts,B,aptr0)
            else error("$T not supported"); end
        else
            if T <: Float32;    ccall((:xfill_32,libknet8),Void,(Cint,Cint,Cfloat, Ptr{Cfloat}, Cint),nrows,ncols,B,aptr0,astep)
            elseif T<: Float64; ccall((:xfill_64,libknet8),Void,(Cint,Cint,Cdouble,Ptr{Cdouble},Cint),nrows,ncols,B,aptr0,astep)
            else error("$T not supported"); end
        end
    else
        length(B) == nelts || throw(DimensionMismatch())
        B = convert(KnetArray{T},B)
        if ncols == 1
            @cuda(cudart,cudaMemcpyAsync,(Ptr{Void},Ptr{Void},Csize_t,UInt32,Ptr{Void}),
                  aptr0, B, nelts*sizeof(T), cudadir(A,B), C_NULL)
        else
            nrows *= sizeof(T); astep *= sizeof(T)
            ccall((:xcopy,libknet8),Void,(Cint,Cint,Ptr{Void},Cint,Ptr{Void},Cint), nrows, ncols, B, nrows, aptr0, astep)
        end
    end
    return A
end

function indexparams{T,N}(A::KnetArray{T,N}, I::Union{Real, UnitRange, Colon}...)
    N > 2 && error("setindex for ndims > 2 not implemented yet")
    skipped = false
    nrows = nelts = 1
    subs1 = ones(Int,2)
    astep = length(A)
    for i=1:length(I)
        Ii = I[i]
        Ai = size(A,i)
        if isa(Ii, Colon)
            Li = Ai
            subs1[i] = 1
        elseif isa(Ii, Real)
            1 <= Ii <= Ai || throw(DimensionMismatch())
            Li = 1
            subs1[i] = Int(Ii)
        else
            1 <= first(Ii) <= last(Ii) <= Ai || throw(DimensionMismatch())
            Li = length(Ii)
            subs1[i] = first(Ii)
        end
        nelts *= Li
        if !skipped
            nrows *= Li
            if Li < Ai
                skipped = true
                astep = stride(A,i+1)
            end
        end
    end
    ncols = div(nelts, nrows)
    firstindex = sub2ind(size(A),subs1...)
    return (nelts,nrows,ncols,firstindex,astep)
end

# Concatenation:
import Base: hcat, vcat, cat

# Benchmarks in μs for hcat and vcat: a=rand(1000,1000) v=rand(1000), t=v'
#		cpu	gpu	g->c->g	vkernel
# hcat(a,a)	2350	225	16160
# hcat(a,v)	1230	115	6490
# hcat(v,a)	1220	120	6490
# hcat(v,v)	3.53	12.53	48.49
# vcat(a,a)	2630	10980	16590	665
# vcat(a,t)	1350	10860	6550	338
# vcat(t,a)	1360	10850	6570	338
# vcat(v,v)	2.13	12.33	45.40	13.58

# setindex! methods called by hcat/vcat:
# hcat(v,v): I = (Colon(),1:1) I = (Colon(),2:2)
# vcat(v,v): uses single index
# hcat(m,m): I = (Colon(),1:5) I = (Colon(),6:10)
# vcat(m,m): I = (1:3,Colon()) I = (4:6,Colon())

function hcat{T}(a::KnetVecOrMat{T}, b::KnetVecOrMat{T})
    size(a,1)==size(b,1) || throw(DimensionMismatch())
    c1 = size(a,1)
    c2 = size(a,2) + size(b,2)
    c = KnetArray(T, (c1,c2))
    c[:,1:size(a,2)] = a
    c[:,1+size(a,2):end] = b
    return c
end

function vcat{T}(a::KnetVector{T}, b::KnetVector{T})
    c = KnetArray(T, length(a)+length(b))
    c[1:length(a)] = a
    c[1+length(a):end] = b
    return c
end

function vcat{T}(a::KnetVecOrMat{T}, b::KnetVecOrMat{T})
    size(a,2)==size(b,2) || throw(DimensionMismatch())
    c1 = size(a,1) + size(b,1)
    c2 = size(a,2)
    c = KnetArray(T, (c1,c2))
    c[1:size(a,1),:] = a
    c[1+size(a,1):end,:] = b
    return c
end

function cat{T}(d, a::KnetVecOrMat{T}, b::KnetVecOrMat{T})
    if     d==1; vcat(a,b)
    elseif d==2; hcat(a,b)
    else error("cat($d) not implemented.")
    end
end

# Utilities:

# Generalizing low level copy using linear indexing to/from gpu arrays:
# copy!{T}(dest::Array{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer)
# Note that this is an unsafe operation, no argument or bounds checking performed.
# Defined in Base:
# unsafe_copy!{T}(dest::Ptr{T}, src::Ptr{T}, n) at array.jl:73
# unsafe_copy!{T}(dest::Array{T,N}, doffs, src::Array{T,N}, soffs, n) at array.jl:79

import Base: unsafe_copy!, copy

function unsafe_copy!{T}(dest::Union{KnetArray{T},Array{T}}, doffs, src::Union{KnetArray{T},Array{T}}, soffs, n; stream=C_NULL)
    @cuda(cudart,cudaMemcpyAsync,(Ptr{Void},Ptr{Void},Csize_t,UInt32,Ptr{Void}),
          pointer(dest,doffs), pointer(src,soffs), n*sizeof(T), cudadir(dest,src), stream)
    return dest
end

function cudadir(a,b)
    deva = isa(a,KnetArray) && a.ptr.dev >= 0
    devb = isa(b,KnetArray) && b.ptr.dev >= 0
    if !deva && !devb; return 0
    elseif deva && !devb; return 1
    elseif !deva && devb; return 2
    elseif deva && devb;  return 3
    end
end

copy(a::KnetArray)=unsafe_copy!(similar(a),1,a,1,length(a))

# Efficient fill:
for S in (32,64); T = Symbol("Float$S"); F = "fill_$S"
    @eval function knetfill!(a::KnetArray{$T},v::$T,off,len)
        ccall(($F,$libknet8),Void,(Cint,$T,Ptr{$T}),len,v,pointer(a,off))
    end
end

# AutoGrad functions:
import AutoGrad: zeroslike, sum_outgrads, OneHot, unary_nd, indexed_function, isequivalent
zeroslike(a::KnetArray)=zeros(a)
sum_outgrads{T}(a::KnetArray{T},b::KnetArray{T})=(a+b)
sum_outgrads(a::KnetArray,b::OneHot)=setindex!(a,sum_outgrads(getindex(a,b.index...),b.value),b.index...)
unary_nd(f, x::KnetArray, eps) = reshape(eltype(x)[unary_nd(indexed_function(f, x, i), x[i], eps) for i in 1:length(x)], size(x))
isequivalent(x::Union{KnetArray,AbstractArray}, y::Union{KnetArray,AbstractArray}; o...)=(length(x)==length(y) && all(i->isequivalent(x[i],y[i];o...), 1:length(x)))

# Hack for printing without copying the whole KnetArray and without inheriting AbstractArray:
import Base: display, summary
type KnetDisplay{T,N} <: AbstractArray{T,N}; a::KnetArray{T,N}; end
getindex(a::KnetDisplay, i...) = getindex(a.a, i...)
size(a::KnetDisplay) = size(a.a)
summary(a::KnetDisplay) = summary(a.a)
summary(a::KnetArray) = string(Base.dims2string(size(a)), " ", typeof(a))
display(a::KnetArray) = display(KnetDisplay(a))


text:
@testset "Scenario" begin

    # Test Scenario
    T = Float64
    N = 10
    p = 3
    model = UnicycleGame(p=p)
    roadway_opts = HighwayRoadwayOptions()
    roadway = build_roadway(roadway_opts)
    player = Vector{Player}(undef, p)
    player[1] = Player(model, roadway.lane[1])
    player[2] = Player(model, roadway.lane[2])
    player[3] = Player(model, roadway.lane[3])
    sce = Scenario(model, roadway, player)
    @test typeof(sce) <: Scenario

    x0 = zeros(model.n)
    x0[4] = roadway_opts.lane_width/2
    x0[5] = -roadway_opts.lane_width/2
    x0[6] = -roadway_opts.lane_width/2
    @test get_state(sce) == SVector{model.n}(x0)

    @test all(get_control_bounds(sce) .==
        (-Inf*ones(SVector{model.m,T}), Inf*ones(SVector{model.m,T})))

    game_obj = GameObjective(N, sce)
    @test typeof(game_obj) <: Algames.GameObjective

    game_con = GameConstraintValues(N, sce)
    @test typeof(game_con) <: Algames.GameConstraintValues

    dt = 0.1
    opts = Options()
    prob = GameProblem(N, dt, sce, opts)
    @test typeof(prob) <: Algames.GameProblem

end


text:
mutable struct polynomial_matrix_t <: LCMType
    timestamp::Int64
    rows::Int32
    cols::Int32
    polynomials::Matrix{polynomial_t}
end

@lcmtypesetup(polynomial_matrix_t,
    polynomials => (rows, cols),
)


text:
using Flux
using Zygote
using MAT
#using Plots
using Statistics
using Distributions
using LinearAlgebra
using Random
using BSON
using JLD
using PyPlot
using CUDA
#using Base.Iterators: partition

gpu_id = 1
CUDA.device!(gpu_id)

Xtrain = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","Xtrain")
Ytrain = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","Ytrain")
XValid = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","XValid")
YValid = JLD.load("Data/XYtrain_XYvalid_processed_T1w.jld","YValid")

idx_shuffle = shuffle(collect(1:size(Xtrain,4)))
Xtrain = Xtrain[:,:,:,idx_shuffle]
Ytrain = Ytrain[:,:,:,idx_shuffle]

Xtrain = reshape(Xtrain,size(Xtrain,1),size(Xtrain,1),2,size(Xtrain,4))

Xtrain = gpu(Xtrain)
Ytrain = gpu(Ytrain)

batch_size = 10

train_loader = Flux.Data.DataLoader(Xtrain,Ytrain, batchsize=batch_size,shuffle=true)

#XValid = XValid[:,:,1,:]
XValid = reshape(XValid,size(XValid,1),size(XValid,1),2,size(XValid,4))

XValid = gpu(XValid)
YValid = gpu(YValid)

channels_in=2
channels_out=1

# Unet define
function _random_normal(shape...)
  return Float32.(rand(Normal(0.f0,0.02f0),shape...))
end

BatchNormWrap(out_ch) = BatchNorm(out_ch)

ResNetConvBlock(chs, kernel = (3, 3)) =
    Chain(Conv(kernel, chs=>chs,pad = (1, 1);init=_random_normal),
          BatchNorm(chs),
          x->relu.(x),
          Conv(kernel, chs=>chs,pad = (1, 1);init=_random_normal),
          BatchNorm(chs))

r = Chain(
            Conv((4,4),channels_in=>64,pad=(1,1),stride=(2,2);init=_random_normal),
            x->relu.(x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            SkipConnection(ResNetConvBlock(64),(mx,x)->mx+x),
            ConvTranspose((2,2),64=>channels_out,stride=(2,2);init=_random_normal),
            BatchNorm(channels_out)
        )

r = gpu(r)

#loss(x,y) = mean(abs,(u(x).-y))
loss(x,y) = mean((r(x).-y).^2)
validationloss() = mean(Flux.mse(r(XValid[:,:,:,(idx-1)*100+1:idx*100]),YValid[:,:,:,(idx-1)*100+1:idx*100]) for idx=1:round(Int,size(XValid,4)/100))
#validationloss() = mean(mean(abs,(u(XValid[:,:,:,(idx-1)*100+1:idx*100]).-YValid[:,:,:,(idx-1)*100+1:idx*100])) for idx=1:round(Int,size(XValid,4)/100))

opt = ADAM(0.001)

function my_custom_train!(loss, ps, data, opt)
    # training_loss is declared local so it will be available for logging outside the gradient calculation.
    local training_loss
    ps = Params(ps)
    for d in data
      if d isa AbstractArray{<:Number}
        gs = Flux.gradient(ps) do
          training_loss = loss(d)
          return training_loss
        end
      else
        gs = Flux.gradient(ps) do
          training_loss = loss(d...)
          return training_loss
        end
      end
      # Insert whatever code you want here that needs training_loss, e.g. logging.
      # logging_callback(training_loss)
      # Insert what ever code you want here that needs gradient.
      # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.
      #println("training_loss = ",training_loss,"; gradient[1] = ",gs[ps[1]],"; W1 = ",ps.order.data[1][1],"; b1 = ",ps.order.data[2][1])
      println("training_loss = ",training_loss)
      Flux.update!(opt, ps, gs)
      # Here you might like to check validation set accuracy, and break out to do early stopping.
    end
  end

loss_valid = zeros(100)
for epoch in 1:100
  learn_rate = 0.001*exp(log(0.1)/99*(epoch-1))
  opt = ADAM(learn_rate)
  #@time Flux.train!(loss, Flux.params(u), train_loader, opt, cb = cb)
  @time my_custom_train!(loss, Flux.params(r), train_loader, opt)
  #@time Flux.train!(loss, Flux.params(u), train_set, opt)
  #model = cpu(u)
  #BSON.@save string("Results/model_resnet_ptx_processed_7dir_gpu_epoch",string(epoch),".bson") model
  #u = gpu(model)
  loss_valid[epoch] = validationloss()
end

model = cpu(r)
BSON.@save string("Results/model_resnet_ptx_processed_7Vol_gpu_T1w_epoch100.bson") model

#=
loss_valid = zeros(50)
for epoch in 1:50
  learn_rate = 0.001*exp(log(0.1)/99*(epoch-1))
  opt = ADAM(learn_rate)
  @time Flux.train!(loss, Flux.params(r), train_loader, opt)
  #@time my_custom_train!(loss, Flux.params(r), train_loader, opt)
  #@time Flux.train!(loss, Flux.params(u), train_set, opt)
  #model = cpu(u)
  #BSON.@save string("Results/model_resnet_ptx_processed_7dir_gpu_epoch",string(epoch),".bson") model
  #u = gpu(model)
  loss_valid[epoch] = validationloss()
end

model = cpu(r)
BSON.@save string("Results/model_resnet_ptx_processed_7Vol_gpu_T1w_epoch50.bson") model
=#

#test
#=
XValid = JLD.load("XYtrain_XYvalid_processed.jld","XValid")
YValid = JLD.load("XYtrain_XYvalid_processed.jld","YValid")
=#

BSON.@load string("Results/model_resnet_ptx_processed_7Vol_gpu_T1w_epoch50.bson") model
slc_list = collect(1:100)
dir_list = collect(1:7)
XtestOutput4write = zeros(size(XValid,1),size(XValid,2),length(slc_list),length(dir_list))

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),2,1)
Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
for idx_dir_test in 1:length(dir_list)
  for idx_slc_test in 1:length(slc_list)
    idx_test = idx_slc_test + (idx_dir_test-1)*100
    Xtest[:,:,:,1] = XValid[:,:,:,idx_test]
    #Ytest[:,:,1,1] = YValid[:,:,idx_slc_test,idx_dir_test]
    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = model(Xtest)
  end
end

file = matopen("Results/XValidOutput_subj5_processed_7Vol_gpu_resnet_T1w_epoch50.mat", "w")
write(file, "XtestOutput4write", XtestOutput4write)
close(file)

# testing on training data
BSON.@load string("model_unet_ptx_processed_7dir_epoch",string(10),".bson") u
XValid = Xtrain[:,:,:,1:700]
XValid = reshape(XValid,size(XValid,1),size(XValid,2),100,7)
XtestOutput4write = zeros(size(XValid))

Xtest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
Ytest = zeros(Float32,size(XValid,1),size(XValid,2),1,1)
for idx_dir_test in 1:size(XValid,4)
  for idx_slc_test in 1:size(XValid,3)
    Xtest[:,:,1,1] = XValid[:,:,idx_slc_test,idx_dir_test]
    #Ytest[:,:,1,1] = YValid[:,:,idx_slc_test,idx_dir_test]
    XtestOutput4write[:,:,idx_slc_test,idx_dir_test] = u(Xtest)
  end
end

file = matopen("XValidOutput_traingsubj1_processed_7dir_epoch10.mat", "w")
write(file, "XtestOutput4write", XtestOutput4write)
close(file)


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

# directory of CBLIB files
cblib_dir = joinpath(ENV["HOME"], "cblib/cblib.zib.de/download/all")
if !isdir(cblib_dir)
    @warn("CBLIB download folder not found")
    cblib_dir = joinpath(@__DIR__, "cblib_data")
    cblib_diverse = String[]
else
    cblib_diverse = [
        "expdesign_D_8_4", # psd, exp
        "port_12_9_3_a_1", # psd, soc, exp
        "tls4", # soc
        "ck_n25_m10_o1_1", # rsoc
        "rsyn0805h", # exp
        "2x3_3bars", # psd
        "HMCR-n20-m400", # power
        "classical_20_0", # soc, orthant
        "achtziger_stolpe06-6.1flowc", # rsoc, orthant
        "LogExpCR-n100-m400", # exp, orthant
        ]
end

relaxed_tols = (default_tol_relax = 1000,)
insts = OrderedDict()
insts["minimal"] = [(("expdesign_D_8_4",), nothing, relaxed_tols)]
insts["fast"] = [((inst,), nothing, relaxed_tols) for inst in cblib_diverse]
insts["various"] = insts["fast"]
return (CBLIBJuMP, insts)

# other sets:

# cblib_power_small = [
#     "HMCR-n20-m400",
#     "HMCR-n100-m400",
#     "HMCR-n500-m400",
#     "HMCR-n20-m800",
#     "HMCR-n100-m800",
#     "HMCR-n500-m800",
#     ]
#
# cblib_exp_small = [
#     "bss1",
#     "demb782",
#     "bss2",
#     "demb781",
#     "synthes1",
#     "gptest",
#     "rijc781",
#     "synthes2",
#     "rijc784",
#     "rijc785",
#     "rijc786",
#     "syn05m",
#     "rijc782",
#     "batchdes",
#     "synthes3",
#     "rijc783",
#     "syn10m",
#     "syn05h",
#     "beck751",
#     "beck752",
#     "beck753",
#     "fiac81b",
#     "batch",
#     "syn05m02m",
#     "syn15m",
#     "fang88",
#     "syn10h",
#     "demb761",
#     "demb762",
#     "demb763",
#     "syn20m",
#     "syn05m02h",
#     "syn05m03m",
#     "fiac81a",
#     "syn30m",
#     "syn10m02m",
#     "rijc787",
#     "syn15h",
#     "ravem",
#     "syn05m04m",
#     "enpro56",
#     "syn05m03h",
#     "syn20h",
#     "syn40m",
#     "rsyn0805m",
#     "enpro48",
#     "syn10m02h",
#     "syn15m02m",
#     "syn10m03m",
#     "rsyn0810m",
#     "syn05m04h",
#     "rsyn0815m",
#     "syn30h",
#     "rsyn0820m",
#     "syn20m02m",
#     "rsyn0805h",
#     "syn10m04m",
#     "rsyn0830m",
#     "isil01",
#     "syn10m03h",
#     "syn15m02h",
#     "syn15m03m",
#     "rsyn0810h",
#     "syn40h",
#     "rsyn0840m",
#     "rsyn0815h",
#     "car",
#     "syn30m02m",
#     "syn20m02h",
#     "rsyn0820h",
#     "syn20m03m",
#     "syn10m04h",
#     "rsyn0805m02m",
#     "syn15m04m",
#     "rsyn0830h",
#     "rsyn0810m02m",
#     "syn15m03h",
#     "batchs101006m",
#     "syn40m02m",
#     "gp_dave_1",
#     "rsyn0815m02m",
#     "rsyn0840h",
#     "syn20m04m",
#     "syn30m02h",
#     "syn30m03m",
#     "rsyn0820m02m",
#     "syn20m03h",
#     "rsyn0805m03m",
#     "rsyn0805m02h",
#     "syn15m04h",
#     "jha88",
#     "batchs121208m",
#     "rsyn0830m02m",
#     "rsyn0810m02h",
#     "rsyn0810m03m",
#     "syn40m02h",
#     "syn40m03m",
#     "batchs151208m",
#     "syn30m04m",
#     "syn20m04h",
#     "varun",
#     "rsyn0840m02m",
#     "rsyn0815m02h",
#     ]
#
# cblib_exp_medium = [
#     "rsyn0815m03m",
#     "syn30m03h",
#     "rsyn0805m04m",
#     "gp_dave_2",
#     "rsyn0820m03m",
#     "rsyn0820m02h",
#     "rsyn0805m03h",
#     "batchs201210m",
#     "rsyn0810m04m",
#     "syn40m04m",
#     "rsyn0830m02h",
#     "rsyn0830m03m",
#     "rsyn0810m03h",
#     "rsyn0815m04m",
#     "syn40m03h",
#     "syn30m04h",
#     "LogExpCR-n20-m400",
#     "gp_dave_3",
#     "rsyn0820m04m",
#     "rsyn0815m03h",
#     "rsyn0840m02h",
#     "rsyn0840m03m",
#     "LogExpCR-n100-m400",
#     "rsyn0805m04h",
#     "rsyn0820m03h",
#     "LogExpCR-n500-m400",
#     "rsyn0810m04h",
#     "rsyn0830m04m",
#     "syn40m04h",
#     "rsyn0830m03h",
#     "rsyn0815m04h",
#     "rsyn0840m04m",
#     "rsyn0820m04h",
#     "rsyn0840m03h",
#     "mra01",
#     "rsyn0830m04h",
#     ]


text:
module SelectivelyAdaptiveLasso

# temp for dev
export Bases, basis_search
export coordinate_descent
# eventual exports
export SALSpec, SALFit, fit, predict

include("bases.jl")
include("ccd.jl")

import StatsBase: sample
import Statistics: mean

mutable struct SALSpec
    # stopping hyperparameters
    λ::Float64 # desired regularization strength
    max_iter::Int # maximum iterations to step through

    # fitting hyperparameters
    bases_per_iter::Int # number of bases to add at once
    subsample_pct::Float64
    subsample_n::Union{Int, Nothing}
    feat_pct::Float64
    tol::Float64 # tolerance for lasso coordinate descent convergence (loss = mse + penalty)
end

function SALSpec(;
        λ::Real = 0.01,
        max_iter::Int = 5000,
        bases_per_iter::Int = 1,
        subsample_pct::Float64 = 0.1,
        subsample_n::Union{Int, Nothing} = nothing,
        feat_pct::Float64 = 1,
        tol::Real = 1e-4,
    )
    return SALSpec(
        float(λ), max_iter, bases_per_iter,
        subsample_pct, subsample_n, feat_pct,
        float(tol),
    )
end

mutable struct SALFit
    β::Dict{BasisIndex, Float64}
end

function predict(sal_fit::SALFit, X::Matrix{Float64})
    bases = Bases(
        Features([FeatureVector(x) for x in eachcol(X)]),
        indices=keys(sal_fit.β)
    )
    return bases * sal_fit.β
end

function fit(
    # model object and hyperparameters
        sal::SALSpec,
    # training data
        X::Matrix{Float64},
        Y::Vector{Float64};
    # validation data for early stopping
        X_val::Union{Nothing, Matrix{Float64}} = nothing,
        Y_val::Union{Nothing, Vector{Float64}} = nothing,
    # logging
        verbose::Bool = false,
        print_iter::Int = 100,
    )::Tuple{SALFit, Any}
    #=
    Fits the selectively adaptive lasso.
    =#

    val = (!isnothing(X_val) & !isnothing(Y_val))
    n,p = size(X)
    X = [FeatureVector(x) for x in eachcol(X)]
    λ = sal.λ * max(sum(Y[Y.>0]), sum(Y[Y.<0])) # makes it so λ=1 => β=0
    ρ = 0

	bases = Bases(X)
	β, R, _ = coordinate_descent(bases, Y, λ=λ, tol=sal.tol)

    mse = [mean(R.^2)]
    if val
        X_val = [FeatureVector(x) for x in eachcol(X_val)]
        bases_val = Bases(X_val)
        R_val = Y_val - bases_val*β
        mse_val = [mean(R_val.^2)]
    end

    if isnothing(sal.subsample_n)
        subsample_n = min(Int(ceil(sal.subsample_pct*n)), n)
    else
        subsample_n = sal.subsample_n
    end
    feat_n = min(Int(ceil(sal.feat_pct*p)), p)

    for i in 1:sal.max_iter

    	for j in 1:sal.bases_per_iter
            index, basis = basis_search(
                X, R, λ,
                subsample_n=subsample_n, feat_n=feat_n,
                off_limits = deepcopy(keys(bases))
            )
	        add_basis!(bases, index, basis)
            if val
                basis_val = build_basis(X_val, index)
                add_basis!(bases_val, index, basis_val)
            end
	    end

        β, R, l = coordinate_descent(bases, Y, λ=λ, β=β, tol=sal.tol)
        filter_bases!(bases, keys(β))

        push!(mse, mean(R.^2))
        if val
            R_val = Y_val - bases_val*β
            push!(mse_val, mean(R_val.^2))
        end

        if verbose & (i % print_iter == 0)
		    print((sqrt(mse[end]), sqrt(mse_val[end]), length(β), length(bases.dict)))
		    print("\n")
		end
    end

    return SALFit(β), (mse, mse_val)
end

end # module


text:
# ---
# title: 1292. Maximum Side Length of a Square with Sum Less than or Equal to Threshold
# id: problem1292
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/maximum-side-length-of-a-square-with-sum-less-than-or-equal-to-threshold/description/>
# hidden: true
# ---
#
# Given a `m x n` matrix `mat` and an integer `threshold`. Return the maximum
# side-length of a square with a sum less than or equal to `threshold` or return
# **0** if there is no such square.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/05/e1.png)
#
#
#
#     Input: mat = [[1,1,3,2,4,3,2],[1,1,3,2,4,3,2],[1,1,3,2,4,3,2]], threshold = 4
#     Output: 2
#     Explanation: The maximum side length of square with sum less than 4 is 2 as shown.
#
#
# **Example 2:**
#
#
#
#     Input: mat = [[2,2,2,2,2],[2,2,2,2,2],[2,2,2,2,2],[2,2,2,2,2],[2,2,2,2,2]], threshold = 1
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: mat = [[1,1,1,1],[1,0,0,0],[1,0,0,0],[1,0,0,0]], threshold = 6
#     Output: 3
#
#
# **Example 4:**
#
#
#
#     Input: mat = [[18,70],[61,1],[25,85],[14,40],[11,96],[97,96],[63,45]], threshold = 40184
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= m, n <= 300`
#   * `m == mat.length`
#   * `n == mat[i].length`
#   * `0 <= mat[i][j] <= 10000`
#   * `0 <= threshold <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module m1testsingle4
using FinEtools
using FinEtoolsDeforNonlinear
using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
using FinEtoolsDeforNonlinear.FEMMDeforNonlinearBaseModule: stiffness
using LinearAlgebra: norm, cross
using SparseArrays
using DelimitedFiles
# using Debugger
using Interpolations
using Test
function test()
    mr = DeforModelRed3D
    E, nu = 7.0*phun("MPa"), 0.3
    m = MatDeforNeohookean(mr, E, nu)
    L= 6/2*phun("mm");
    H = 2/2*phun("mm");
    W = 2/2*phun("mm");
    umag=-1.5*phun("mm");# Magnitude of the displacement
    nincr = 48

    fens, fes = H8hexahedron([0 0 0; 1.2 0 -0.1; 1.0 0.9 0.1; -0.05 1.1 0; 0 0 1.03; 1.2 0 1.1; 1.0 0.95 0.81; -0.05 1.01 0.9], 2, 1, 1)
    femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

    geom = NodalField(fens.xyz)
    un1 = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field
    un = NodalField(zeros(size(fens.xyz, 1), 3)) # displacement field

    # Clamped cross-section
    node_list = selectnode(fens, box = [0,0,-Inf,Inf,-Inf,Inf], inflate = W/1000);
    setebc!(un, node_list, true, 1, 0.0)
    node_list = selectnode(fens, box = [-Inf,Inf,0,0,-Inf,Inf], inflate = W/1000);
    setebc!(un, node_list, true, 2, 0.0)
    node_list = selectnode(fens, box = [-Inf,Inf,-Inf,Inf,0,0], inflate = W/1000);
    setebc!(un, node_list, true, 3, 0.0)
    table = LinearInterpolation([0, 0.25, 0.5, 0.75,1], umag*[0,-1,0,1,0])
    move(lambda) = table(lambda);
    movingl = selectnode(fens, box = [L,L,-Inf,Inf,-Inf,Inf], inflate = W/1000);
    setebc!(un, movingl, true, 1, 0.0)
    applyebc!(un)
    numberdofs!(un)
    copyto!(un1, un);

    tn, dtn = 0.0, 1.0
    load_multipliers = (0:1:nincr)/nincr*1.0;
    load_increments = diff(unique(sort(load_multipliers)));

    # For all  load increments
    for incr =1:length(load_increments) # Load-implementation loop

        lambda = sum(load_increments[1:incr];
        dlambda = load_increments[incr];

        # Initial value of the displacement  in the current step
        copyto!(un1, un);

        # Process load-factor-dependent essential boundary conditions:
        setebc!(un1, movingl, true, 1, move(lambda))
        applyebc!(un1)

        # Initialization
        applyebc!(un1) # Apply EBC

        du = deepcopy(un1); # Displacement increment
        du.values[:] = 0.0

        # Initialize the load vector
        F = fill(0.0, un1.nfreedofsyou);

    # If any boundary conditions are inhomogeneous, calculate  the force
    # vector due to the displacement increment. Then update the guess of
    # the new displacement.
    if (aany_nonzero_EBC)
        for i=1:length(model_data.region)
            F = F + nz_ebc_loads(model_data.region{i}.femm, sysvec_assembler, geom, un, unm1, un1-un, dlambda);
        end
        # Provided we got converged results  in the last step, we already
        # have a usable stiffness matrix.
        if (isempty(K)) # we don't have a stiffness matrix
            K=  sparse(un1.nfreedofs,un1.nfreedofs);
            for i=1:length(model_data.region)
                K = K + stiffness(model_data.region{i}.femm, model_data.region{i}.MA, geom, un, unm1, dlambda);
                K = K + stiffness_geo(model_data.region{i}.femm, model_data.region{i}.MA, geom, un, unm1, dlambda);
            end
        end
        un1 = un1 + scatter_sysvec(du, K\F);
    end

    # Iteration loop
    iter=1;
    while true #  Iteration loop

        # Initialize the load vector
        F =0*F;

        # Construct the system stiffness matrix.
        K=  sparse(un1.nfreedofs,un1.nfreedofs);
        for i=1:length(model_data.region)
            K = K + stiffness(model_data.region{i}.femm, model_data.region{i}.MA, geom, un1, un, dlambda);
            K = K + stiffness_geo(model_data.region{i}.femm, model_data.region{i}.MA, geom, un1, un, dlambda);
        end

        # Process the body load
        if (isfield(model_data, 'body_load' ))
            for j=1:length(model_data.body_load)
                body_load =model_data.body_load{j};
                femm = femm_deformation (struct ('material',[],...
                    'fes',body_load.fes,...
                    'integration_rule',body_load.integration_rule));
                fi= force_intensity(struct('magn',body_load.force));
                F = F + distrib_loads(femm, sysvec_assembler, geom, un1, fi, 3);
            end
            clear body_load fi  femm
        end

        # Process the traction boundary condition
        if (isfield(model_data.boundary_conditions, 'traction' ))
            for j=1:length(model_data.boundary_conditions.traction)
                traction =model_data.boundary_conditions.traction{j};
                femm = femm_deformation (struct ('material',[],...
                    'fes',traction.fes,...
                    'integration_rule',traction.integration_rule));
                fi= force_intensity(struct('magn',traction.traction));
                F = F + distrib_loads(femm, sysvec_assembler, geom, un1, fi, 2);
            end
            clear traction fi  femm
        end

        # Process the nodal force boundary condition
        if (isfield(model_data.boundary_conditions, 'nodal_force' ))
            for j=1:length(model_data.boundary_conditions.nodal_force)
                nodal_force =model_data.boundary_conditions.nodal_force{j};
                femm = femm_deformation (struct ('material',[],...
                    'fes',fe_set_P1(struct('conn',reshape(nodal_force.node_list,[],1))),...
                    'integration_rule',point_rule));
                fi= force_intensity(struct('magn',nodal_force.force));
                F = F + distrib_loads(femm, sysvec_assembler, geom, un1, fi, 0);
            end
            clear nodal_force fi femm
        end

        #  External loads vector
        FL = lambda * F;

        # The restoring  force vector
        FR =zeros(un1.nfreedofs,1);
        for i=1:length(model_data.region)
            FR = FR + restoring_force(model_data.region{i}.femm, sysvec_assembler, geom, un1, un, dlambda);
        end

        # Solve the system of linear algebraic equations
        F = (FL + FR);
        dusol= K\F;

        #  Distribute the solution
        du = scatter_sysvec(du, dusol);

        # Do we need line search?
        eta= 1.0;
        if (line_search)
            R0 = dot(F,dusol);
            for  linesrch   =1:1 # How many  searches should we do?
                FR =zeros(un1.nfreedofs,1);
                for i=1:length(model_data.region)
                    FR = FR + restoring_force(model_data.region{i}.femm, sysvec_assembler, geom, un1+eta*du, un, dlambda);
                end
                F = FL + FR;
                R1 = dot(F,dusol);
                a = R0/R1;
                if ( a<0 )
                    eta = a/2 +sqrt((a/2)^2 -a);
                else
                    eta =a/2;
                end
                eta=min( [eta, 1.0] );
            end
        end

        # Increment the displacements
        un1 = un1 + eta*du;

        # Compute the increment data
        maxdu  =eta*max(max(abs(du.values)));
        maxbal=max(abs(F));

        # Report  iteration results?
        if (~isempty(iteration_observer))
            model_data.un = un;
            model_data.un1 = un1;
            model_data.dt = dlambda;
            model_data.maxdu = maxdu;
            model_data.maxbal = maxbal;
            iteration_observer (lambda,iter,du,model_data);
        end

        # Converged?
        if (maxdu <= maxdu_tol)  || (maxbal <= maxbal_tol)
            break; # Breakout of the iteration loop
        end;

        iter=iter+1;

    end #  Iteration loop

    # Update the model data
    model_data.un = un;
    model_data.un1 = un1;
    model_data.dt = dlambda;
    #  Now the reactions
    un1Allfree.values=un1.values;# This field holds the current converged displacements;
    unAllfree=un1Allfree;
    unAllfree.values=un.values;# This field holds the current converged displacements;
    # since all the degrees of freedom are free, the restoring forces
    # can be used to compute the reactions. Note the negative sign: the
    # reactions are the opposite of the resisting forces of the
    # material.
    FR =zeros(un1Allfree.nfreedofs,1);
    for i=1:length(model_data.region)# Compute the resisting forces of the material
        FR = FR - restoring_force(model_data.region{i}.femm, sysvec_assembler, geom, un1Allfree, unAllfree, dlambda);
    end
    model_data.reactions=scatter_sysvec(un1Allfree,FR);

    # Converged
    # Update the FEMMs
    for i=1:length(model_data.region)
       [~,model_data.region{i}.femm]  =...
           restoring_force(model_data.region{i}.femm, sysvec_assembler, ...
           geom, un1, un, dlambda);
    end

    # Report results
   if ~isempty(increment_observer)# report the progress
        increment_observer (lambda,model_data);
    end

    # Reset  the displacement field for the next load step
    unm1 =un;
    un = un1;
    F =0*F;

end # Load-incrementation loop

end

    # K = stiffness(femm, geom, un1, un, tn, dtn)

    # ijs = readdlm("m1testop4.csv", ',', Float64)
    # matK  = sparse(Int.(ijs[:, 1]), Int.(ijs[:, 2]), ijs[:, 3], size(K, 1), size(K, 2))
    # # @show norm(K - matK) / maximum(matK[:])
    # @test norm(K - matK) / maximum(matK[:]) < 1.0e-4
end
end
using .m1testsingle4
m1testsingle4.test()


text:
using LinearAlgebra
import LinearAlgebra.Eigen
using SparseArrays
using BenchmarkTools
using Test
using Distributions
using ELM
R = rand(Uniform(1., 100.), 10000,1000)
MM = [1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0; 1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0; 1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0; 1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0;1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0;1 0 0 0 2; 0 0 3 0 0; 0 0 0 0 0; 0 2 0 0 0 ]

C = R
println("size of matrix: ", size(C) )

SVD =  svd(C)
mySvd =  SVDCM(C)

println("Julia SVD time:")
SVD = @time svd(C)

println("SVD time:")
mySvd = @time SVDCM(C)

F = mySvd.U * mySvd.S * mySvd.V'
#display(F)
@testset "test predictions" begin
    @test F ≈ C atol=1e-5
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    SimplePath(domain)

A simple (or default) path on a spatial `domain`.
"""
struct SimplePath{D<:AbstractDomain} <: AbstractPath{D}
  domain::D
end

Base.iterate(p::SimplePath, state=1) = state > npoints(p.domain) ? nothing : (state, state + 1)


text:
module OpenQuantumSystems

using Reexport
@reexport using QuantumOpticsBase:
    Basis,
    GenericBasis,
    basis,
    tensor,
    ⊗,
    permutesystems,
    @samebases,

    # states.jl
    StateVector,
    Bra,
    Ket,
    basisstate,
    norm,
    dagger,
    normalize,
    normalize!,

    # operators.jl
    AbstractOperator,
    DataOperator,
    expect,
    variance,
    identityoperator,
    ptrace,
    embed,
    dense,
    tr,
    sparse,

    # operators_dense.jl
    Operator,
    DenseOperator,
    DenseOpType,
    projector,
    dm,

    # operators_sparse.jl
    SparseOperator,
    diagonaloperator,
    SparseOpType,

    # superoperators.jl
    AbstractSuperOperator,
    SuperOperator,
    DenseSuperOperator,
    DenseSuperOpType,
    SparseSuperOperator,
    SparseSuperOpType,
    spre,
    spost,

    # fock.jl
    FockBasis,
    number,
    destroy,
    create,
    fockstate,
    coherentstate,
    coherentstate!,
    displace,

    # spin.jl
    SpinBasis,
    sigmax,
    sigmay,
    sigmaz,
    sigmap,
    sigmam,
    spinup,
    spindown,

    # state_definitions.jl
    randstate,
    randoperator,
    thermalstate,
    coherentthermalstate,
    phase_average,
    passive_state,

    # nlevel.jl
    NLevelBasis,
    transition,
    nlevelstate,

    # metrics.jl
    tracenorm,
    tracenorm_h,
    tracenorm_nh,
    tracedistance,
    tracedistance_h,
    tracedistance_nh,
    entropy_vn,
    entropy_renyi,
    fidelity,
    ptranspose,
    PPT,
    negativity,
    logarithmic_negativity,
    entanglement_entropy,

    # pauli.jl
    PauliBasis,
    PauliTransferMatrix,
    DensePauliTransferMatrix,
    ChiMatrix,
    DenseChiMatrix,
    avg_gate_fidelity

using SparseArrays, LinearAlgebra

export
    # operators_dense.jl
    AnnihilationOperator,
    CreationOperator,
    PositionOperator,
    MomentumOperator,
    ShiftOperator,
    OneDenseOperator,

    # superoperators.jl
    Commutator,

    # metrics.jl
    tracenorm,
    tracenorm_h,
    tracenorm_nh,
    tracedistance,
    tracedistance_h,
    tracedistance_nh,

    # molecules.jl
    Mode,
    convert_units,
    franckCondonFactors,
    vibrationalIndices,
    electronicIndices,
    Molecule,
    updateMolecule!,
    updateMolecule,
    getMolStateEnergy,
    getMolStateFC,
    getMolShifts,
    getMolFrequencies,

    # aggregateCore.jl
    AggregateCore,
    getNvib,
    getShifts,
    getFrequencies,
    getAggStateEnergy,

    # aggregateTools.jl
    AggregateTools,
    vibrationalIndices,
    electronicIndices,
    getIndices,
    getIndicesMap,
    getFranckCondonFactors,
    # getFranckCondonFactorsSparse,
    getFCproduct,
    take_el_part,

    # aggregateOperators.jl
    getAggHamSystemSmall,
    getAggHamSystemBig,
    getAggHamBathSmall,
    getAggHamBathBig,
    getAggHamSystemBath,
    getAggHamInteraction,
    getAggHamiltonian,
    # getAggHamiltonianSparse,
    AggregateOperators,

    # aggregate.jl
    Aggregate,
    setupAggregate,
    setupAggregate!,

    # evolution.jl
    get_tspan,
    evolutionOperator,
    evolutionOperatorA,
    evolutionSuperOperator,
    evolutionOperatorArray,
    evolutionSuperOperatorArray,
    evolutionOperatorIterator,
    evolutionSuperOperatorIterator,
    evolutionExact,
    evolutionExact!,
    evolutionApproximate,
    evolutionApproximate!,
    evolution_exact,
    evolution_approximate,
    evolutionOperatorExp,
    evolution_el_part,
    Evolution_SI_exact,
    Evolution_sI_exact,
    Evolution_SS_exact,
    Evolution_sS_exact,

    # schrodinger.jl
    schroedinger,
    schroedinger_dynamic,

    # liouville.jl
    LvN_sI,
    LvN_sS,
    LvN_SI,
    LvN_SS,

    # interaction_picture.jl
    getInteractionHamIPicture,
    getInteractionHamIPictureA,

    # master_exact.jl
    QME_SI_exact,
    QME_SS_exact,
    QME_sI_exact,
    QME_sS_exact,

    # trace.jl
    trace_bath,
    trace_bath_slow,
    get_rho_bath,
    trace_bath_part,
    ad,
    correlation_function,

    # initial_state.jl
    thermal_state,
    thermal_state_old,
    thermal_state_composite,
    thermal_state_composite_old,
    ultrafast_laser_excitation,

    # memory_kernel.jl
    MemoryKernel_1_traced,
    MemoryKernel_2_traced,
    MemoryKernel_3_traced,
    MemoryKernel_4_traced,
    MemoryKernel_traced,

    # rate_constant.jl
    M_aabb_W_bath_intp,
    K_aabb_W_bath_intp,
    M_abcd_W_bath_intp,
    K_abcd_W_bath_intp,

    # master_ansatz.jl
    QME_sI_ansatz_test,
    QME_sI_ansatz_const_int,
    QME_sI_ansatz_const_sch,
    QME_sI_ansatz_linear_sch,
    QME_sI_ansatz_linear2_sch,
    QME_sI_ansatz_upart1_sch,
    QME_sI_ansatz_upart1_int,
    QME_sI_ansatz_upart2_sch,
    QME_sI_ansatz_upart2_int,

    # master_iterative
    W_1_bath,
    normalize_bath,
    QME_sI_iterative,

    # postprocessing.jl
    OperatorVector,
    OperatorVectorArray,
    operator_recast,
    interaction_pic_to_schroedinger_pic,
    schroedinger_pic_to_interaction_pic,
    local_st_to_exciton_st,
    exciton_st_to_local_st,
    tspan_cm_to_fs,

    # scoring.jl
    get_rmse_in_time,
    compare_rho,
    compare_rho_in_time

include("core.jl")
include("operators_dense.jl")
include("superoperators.jl")
include("metrics.jl")
include("molecules.jl")
include("aggregateCore.jl")
include("aggregateTools.jl")
include("aggregateOperators.jl")
include("aggregate.jl")
include("timeevolution_base.jl")
include("evolution.jl")
include("schroedinger.jl")
include("liouville.jl")
include("interaction_picture.jl")
include("master_exact.jl")
include("trace.jl")
include("initial_state.jl")
include("memory_kernel.jl")
include("rate_constant.jl")
include("master_ansatz.jl")
include("master_iterative.jl")
include("postprocessing.jl")
include("scoring.jl")

end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
function e009(N)
    return product_of_triple(N)
end

function solution009()
    N = 1000
    answer = e009(N)
    println("The answer to problem 009 is $answer")
end

function product_of_triple(sum_of_values::Integer)
    for a = 1:sum_of_values
        for b = a+1:Int(floor((sum_of_values-a-1)/2))
            c = sum_of_values - a - b
            if a^2 + b^2 == c^2
                return a*b*c
            end
        end
    end
    return -1
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
function getEdgeLength(args...; Tree::Array{SubString{String},1} = [],
	Start::Int = Int[])
	u = Start + 1
	stop_char = [",", ")", ";"]
  edge_length = Float64[]
	while countnz(Tree[u] .== stop_char) == 0
		u = u + 1
		edge_length = parse(join(Tree[(Start + 1):(u - 1)]))
	end
	lengthStart = (edge_length, u + Start)
	lengthStart = Dict(zip(["edge_length", "end"], lengthStart))
	return lengthStart
end


text:
using FrameFun
using BasisFunctions, Test, DomainSets
using BasisFunctions.Test: test_generic_dict_interface

import BasisFunctions.Test: point_outside_domain
using BasisFunctions.Test: generic_test_discrete_measure, generic_test_measure

point_outside_domain(::ExtensionFrame) = 1.

@testset "ExtensionFrame" begin
    d1 = extensionframe(Fourier(10),0.0..0.5)
    test_generic_dict_interface(d1)
    d2 = extensionframe(Fourier(10)^2,(0.0..0.5)^2)
    @test d1 == extensionframe(0.0..0.5, Fourier(10))
    @test d2 isa TensorProductDict
    @test d2 isa ExtensionFrameTensor
    @test basis(d1) == Fourier(10)
    @test basis(d2) == Fourier(10)^2
    @test support(d1) == 0.0..0.5
    @test support(d2) == (0.0..0.5)^2
    @test measure(d1) isa SubWeight
    @test measure(d2) isa ProductWeight
    @test iscompatible(extensiondual(d1, measure(d1)),component(extensiondual(d2, measure(d2)), 1))
end

@testset "SubWeight" begin
    g = PeriodicEquispacedGrid(3,-1,1)^2
    d = UnitDisk()
    sg = subgrid(g,d)
    μ = discretemeasure(sg)
    @test μ isa DiscreteSubWeight
    generic_test_discrete_measure(μ)
    @test supermeasure(μ) == discretemeasure(g)
    @test !(isnormalized(μ))
    @test restrict(supermeasure(μ),d)≈μ

    g = PeriodicEquispacedGrid(3,-1,1)×PeriodicEquispacedGrid(4,-1,1)
    d = UnitInterval()^2
    @assert subgrid(g,d) isa ProductSubGrid
    sg = subgrid(g,d)
    μ = discretemeasure(sg)
    generic_test_discrete_measure(μ )
    @test μ isa DiscreteTensorSubWeight
    @test supermeasure(μ) ≈ discretemeasure(g)
    @test components(μ) == (component(μ,1),component(μ,2))
    @test component(μ,1) == restrict(component(supermeasure(μ),1),UnitInterval())

    m = JacobiWeight(rand(),rand())
    μ = restrict(m, UnitInterval())
    generic_test_measure(μ)
    @test supermeasure(μ) == m
    μ = submeasure(m, UnitInterval())
    generic_test_measure(μ)
    @test supermeasure(μ) == m
end


text:
"""
GeneralizedGamma(a,b,k)

The *Generalized gamma distribution* with shape parameters `a`, scale `b`, and power `k` has probability density
function

```math
f(x; a, b) = \\frac{(k/b^a) x^{a-1} e^{-(x/b)^k}{\\Gamma(a/k)},
\\quad x > 0
```

```julia
GeneralizedGamma(a, b, k)      # GeneralizedGamma distribution with shape a, scale b, power k

params(d)        # Get the parameters, i.e. (a, b, k)
shape(d)         # Get the shape parameter, i.e. a
scale(b)         # Get the scale parameter, i.e. b
```

External links

* [GeneralizedGamma distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_gamma_distribution)

"""

struct GeneralizedGamma <: ContinuousUnivariateDistribution
    a
    b
    k
end

#### Parameters
params(d::GeneralizedGamma) = (d.a, d.b, d.k)
shape(d::GeneralizedGamma) = d.a
scale(d::GeneralizedGamma) = d.b
power(d::GeneralizedGamma) = d.k

#### Evaluation
function gamma_cdf(a, x)
    return gamma_inc(a,x,0)[1]
end

function gamma_ccdf(a, x)
    return gamma_inc(a,x,0)[2]
end

function cdf(d::Gamma, x::Real)
    if x < zero(x)
        return zero(x)
    else
        shp, scl = params(d)
        return gamma_cdf(shp, x/scl)
    end
end

function ccdf(d::Gamma, x::Real)
    if x < zero(x)
        return one(x)
    else
        shp, scl = params(d)
        return gamma_ccdf(shp, x/scl)
    end
end

function cdf(d::GeneralizedGamma, x::Real)
    if x < zero(x)
        return zero(x)
    else
        shp, scl, pwr = params(d)
        return gamma_cdf(shp/pwr, (x/scl)^pwr)
    end
end

function ccdf(d::GeneralizedGamma, x::Real)
    if x < zero(x)
        return one(x)
    else
        shp, scl, pwr = params(d)
        return gamma_ccdf(shp/pwr, (x/scl)^pwr)
    end
end

function pdf(d::GeneralizedGamma, x::Real)
    if x < zero(x)
        return zero(x)
    else
        shp, scl, pwr = params(d)
        return (pwr/(scl^shp))*x^(shp-1)*exp(-(x/scl)^pwr)/gamma(shp/pwr)
    end
end

function logpdf(d::GeneralizedGamma, x::Real)
    if x < zero(x)
        return zero(x)
    else
        shp, scl, pwr = params(d)
        return log(pwr)-shp*log(scl) + (shp-1)*log(x) - (x/scl)^pwr -loggamma(shp/pwr)
    end
end

function eqcdf(d::GeneralizedGamma, x::Real)
    out = zero(x)
    if x >= zero(x)
        shp, scl, pwr = params(d)
        out = (one(x) + x*gamma(shp/pwr)*gamma_ccdf(shp/pwr, (x/scl)^pwr)/(scl*gamma((shp+1)/pwr)) - gamma_ccdf((shp+1)/pwr, (x/scl)^pwr))
        return out
    end
end

function quantile(d::GeneralizedGamma, p)
    shp, scl, pwr = params(d)
    r = gamma_inc_inv(shp/pwr, p, 1-p)
    return scl*(r^inv(pwr))
end

function rand(rng::AbstractRNG, d::GeneralizedGamma)
    p = rand(rng)
    return quantile(d, p)
end

#### Statistics

function mean(d::GeneralizedGamma)
    shp, scl, pwr = params(d)
    return scl*gamma((shp+1)/pwr)/gamma(shp/pwr)
end

#var(d::Gamma) = d.α * d.θ^2

#skewness(d::Gamma) = 2 / sqrt(d.α)

#kurtosis(d::Gamma) = 6 / d.α

# function mode(d::Gamma)
#     (α, θ) = params(d)
#     α >= 1 ? θ * (α - 1) : error("Gamma has no mode when shape < 1")
# end

# function entropy(d::Gamma)
#     (α, θ) = params(d)
#     α + loggamma(α) + (1 - α) * digamma(α) + log(θ)
# end

# mgf(d::Gamma, t::Real) = (1 - t * d.θ)^(-d.α)

# cf(d::Gamma, t::Real) = (1 - im * t * d.θ)^(-d.α)

# function kldivergence(p::Gamma, q::Gamma)
#     # We use the parametrization with the scale θ
#     αp, θp = params(p)
#     αq, θq = params(q)
#     θp_over_θq = θp / θq
#     return (αp - αq) * digamma(αp) - loggamma(αp) + loggamma(αq) -
#         αq * log(θp_over_θq) + αp * (θp_over_θq - 1)
# end

#### Evaluation & Sampling

# @_delegate_statsfuns Gamma gamma α θ

# gradlogpdf(d::Gamma{T}, x::Real) where {T<:Real} =
#     insupport(Gamma, x) ? (d.α - 1) / x - 1 / d.θ : zero(T)

# function rand(rng::AbstractRNG, d::Gamma)
#     if shape(d) < 1.0
#         # TODO: shape(d) = 0.5 : use scaled chisq
#         return rand(rng, GammaIPSampler(d))
#     elseif shape(d) == 1.0
#         return rand(rng, Exponential(d.θ))
#     else
#         return rand(rng, GammaGDSampler(d))
#     end
# end

# function sampler(d::Gamma)
#     if shape(d) < 1.0
#         # TODO: shape(d) = 0.5 : use scaled chisq
#         return GammaIPSampler(d)
#     elseif shape(d) == 1.0
#         return sampler(Exponential(d.θ))
#     else
#         return GammaGDSampler(d)
#     end
# end

#### Fit model

# struct GammaStats <: SufficientStats
#     sx::Float64      # (weighted) sum of x
#     slogx::Float64   # (weighted) sum of log(x)
#     tw::Float64      # total sample weight

#     GammaStats(sx::Real, slogx::Real, tw::Real) = new(sx, slogx, tw)
# end

# function suffstats(::Type{<:Gamma}, x::AbstractArray{T}) where T<:Real
#     sx = zero(T)
#     slogx = zero(T)
#     for xi = x
#         sx += xi
#         slogx += log(xi)
#     end
#     GammaStats(sx, slogx, length(x))
# end

# function suffstats(::Type{<:Gamma}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Real
#     n = length(x)
#     if length(w) != n
#         throw(DimensionMismatch("Inconsistent argument dimensions."))
#     end

#     sx = zero(T)
#     slogx = zero(T)
#     tw = zero(T)
#     for i = 1:n
#         @inbounds xi = x[i]
#         @inbounds wi = w[i]
#         sx += wi * xi
#         slogx += wi * log(xi)
#         tw += wi
#     end
#     GammaStats(sx, slogx, tw)
# end

# function gamma_mle_update(logmx::Float64, mlogx::Float64, a::Float64)
#     ia = 1 / a
#     z = ia + (mlogx - logmx + log(a) - digamma(a)) / (abs2(a) * (ia - trigamma(a)))
#     1 / z
# end

# function fit_mle(::Type{<:Gamma}, ss::GammaStats;
#     alpha0::Float64=NaN, maxiter::Int=1000, tol::Float64=1e-16)

#     mx = ss.sx / ss.tw
#     logmx = log(mx)
#     mlogx = ss.slogx / ss.tw

#     a::Float64 = isnan(alpha0) ? (logmx - mlogx)/2 : alpha0
#     converged = false

#     t = 0
#     while !converged && t < maxiter
#         t += 1
#         a_old = a
#         a = gamma_mle_update(logmx, mlogx, a)
#         converged = abs(a - a_old) <= tol
#     end

#     Gamma(a, mx / a)
# end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


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

@doc raw"""
	load_cap_reserve_margin(setup::Dict, path::AbstractString, inputs_crm::Dict)

Function for reading input parameters related to planning reserve margin constraints
"""
function load_cap_reserve_margin(setup::Dict, path::AbstractString, inputs_crm::Dict)
    # Definition of capacity reserve margin (crm) by locational deliverability area (LDA)
    inputs_crm["dfCapRes"] = DataFrame(
        CSV.File(joinpath(path, "Capacity_reserve_margin.csv"), header = true),
        copycols = true,
    )

    # Ensure float format values:

    # Identifying # of planning reserve margin constraints for the system
    res = count(s -> startswith(String(s), "CapRes"), names(inputs_crm["dfCapRes"]))
    first_col = findall(s -> s == "CapRes_1", names(inputs_crm["dfCapRes"]))[1]
    last_col = findall(s -> s == "CapRes_$res", names(inputs_crm["dfCapRes"]))[1]
    inputs_crm["dfCapRes"] = Matrix{Float64}(inputs_crm["dfCapRes"][:, first_col:last_col])
    inputs_crm["NCapacityReserveMargin"] = res

    println("Capacity_reserve_margin.csv Successfully Read!")

    return inputs_crm
end

@doc raw"""
	load_cap_reserve_margin_trans(setup::Dict, inputs_crm::Dict, network_var::DataFrame)

Function for reading input parameters related to participation of transmission imports/exports in capacity reserve margin constraint.
"""
function load_cap_reserve_margin_trans(setup::Dict, inputs_crm::Dict, network_var::DataFrame)

    res = inputs_crm["NCapacityReserveMargin"]

    first_col_trans = findall(s -> s == "CapRes_1", names(network_var))[1]
    last_col_trans = findall(s -> s == "CapRes_$res", names(network_var))[1]
    dfTransCapRes = network_var[:, first_col_trans:last_col_trans]
    inputs_crm["dfTransCapRes"] = Matrix{Float64}(dfTransCapRes[completecases(dfTransCapRes), :])

    first_col_trans_derate = findall(s -> s == "DerateCapRes_1", names(network_var))[1]
    last_col_trans_derate = findall(s -> s == "DerateCapRes_$res", names(network_var))[1]
    dfDerateTransCapRes = network_var[:, first_col_trans_derate:last_col_trans_derate]
    inputs_crm["dfDerateTransCapRes"] = Matrix{Float64}(dfDerateTransCapRes[completecases(dfDerateTransCapRes), :])

    first_col_trans_excl = findall(s -> s == "CapRes_Excl_1", names(network_var))[1]
    last_col_trans_excl = findall(s -> s == "CapRes_Excl_$res", names(network_var))[1]
    dfTransCapRes_excl = network_var[:, first_col_trans_excl:last_col_trans_excl]
    inputs_crm["dfTransCapRes_excl"] = Matrix{Float64}(dfTransCapRes_excl[completecases(dfTransCapRes_excl), :])

    return inputs_crm
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

"""
    true_obliquity(jd) -> t_eps

### Purpose ###

Return the true obliquity of the ecliptic for a given Julian date

### Explanation ###

The function is used by the [`co_aberration`](@ref) procedure.

### Arguments ###

* `jd`: Julian date.

### Output ###

* `t_eps`: true obliquity of the ecliptic, in radians

### Example ###

```jldoctest
julia> using AstroLib

julia> true_obliquity(jdcnv(1978,01,7,11, 01))
0.4090953896211926
```

### Notes ###

The function calls [`mean_obliquity`](@ref).
"""
function true_obliquity(jd::Real)
    eps0 = mean_obliquity(jd)
    t_eps = eps0 + sec2rad(nutate(jd)[2])
    return t_eps
end


text:
"""
    \\input{qualifier}{rpath}

Resolve a command of the form `\\input{qualifier}{rpath}` where

* `qualifier` indicates what should be inserted (either a code block or a
   generated figure).
* `rpath` indication where it is (in Unix format, see [`parse_rpath`](@ref)).

This command is expected to be used to insert code blocks or the output of
code blocks.

Different actions can be taken based on the qualifier:
1. `{lang}` will insert the code in the script with appropriate language
    indicator. If you want to just include it as plain text, use `{plaintext}`.
3. `{plot}` or `{plot:id}` will look for a displayable image file (gif, png,
    jp(e)g or svg) in `[path_to_script_dir]/output/` and will add an `img`
    block as a result.
    If an `id` is specified, it will try to find an image with the same root
    name ending with `id.ext` where `id` can help identify a specific image if
    several are generated by the script, typically a number will be used,
    e.g.: `\\input{plot:4}{ex2}`.
"""
function lx_input(lxc::LxCom, _)
    qualifier = lowercase(stent(lxc.braces[1]))
    rpath     = stent(lxc.braces[2])
    # check the qualifier
    if startswith(qualifier, "plot")
        # check if id is given e.g. \input{plot:5}{ex2}
        s  = split(qualifier, ":")
        id = length(s) > 1 ? s[2] : ""
        return _lx_input_plot(rpath, id)
    end
    # assume it's a language e.g. \input{julia}{ex2}
    return _lx_input_code(rpath; lang=qualifier)
end

"""Helper function to input a script."""
function _lx_input_code(rpath::AS; lang="")::String
    code = ""
    try
        fp,  = resolve_rpath(rpath, lang)
        code = read(fp, String)
    catch e
        return html_err(e.m)
    end
    return html_code(code, lang)
end

"""Helper function to input a plot generated by a script."""
function _lx_input_plot(rpath::AS, id::AS="")
    cp    = form_codepaths(rpath)
    fname = splitext(cp.script_name)[1]
    odir  = cp.out_dir
    pname = fname * id
    # find a plt in output
    isdir(odir) || return html_err("Couldn't find an output directory " *
                                   "associated with '$rpath' when trying " *
                                   "to input a plot.")
    for (root, _, files) in walkdir(odir)
        for (f, e) ∈ splitext.(files)
            lc_e = lowercase(e)
            if f == pname && lc_e ∈ (".gif", ".jpg", ".jpeg", ".png", ".svg")
                # construct a relative path to the plot
                reldir = odir[length(PATHS[:site])+1:end]
                ppath  = unixify(joinpath(reldir, pname * lc_e))
                return html_img(ppath)
            end
        end
    end
    return html_err("Couldn't find a relevant image when trying to input " *
                    "a plot relative to '$rpath'.")
end

"""
    \\output{...}

Return the output of a code. Possibly including the result of the code
(`res=true`) and possibly reprocessing the whole (`reproc=true`).
At most one will be true. See [`lx_show`](@ref).
"""
function lx_output(lxc::LxCom, lxd::Vector{LxDef};
                   reproc::Bool=false, res::Bool=false)
    rpath   = stent(lxc.braces[1])
    cpaths  = form_codepaths(rpath)
    outpath = cpaths.out_path
    respath = cpaths.res_path
    # does output exist?
    isfile(outpath) || return html_err("`$(lxc.ss)`: could not find the " *
                                       "relevant output file.")
    output = read(outpath, String)
    # if result should be appended, check the relevant file and append
    if res
        isfile(respath) || return html_err("`$(lxc.ss)`: could not find " *
                                           "the relevant result file.")
        result = read(respath, String)
        if result != "nothing"
            if !isempty(output) && !endswith(output, "\n")
                output *= "\n"
            end
            output *= replace(result, r"Main\.FD_SANDBOX_[0-9]+\." => "")
        end
    end
    # should it be reprocessed ?
    reproc || return html_code(output; class="code-output")
    return reprocess(output, lxd)
end

"""
    \\textoutput{...}

Same as [`output`](@ref) but with re-processing.
"""
lx_textoutput(lxc::LxCom, lxd) = lx_output(lxc, lxd; reproc=true)

"""
    \\show{...}

Same as [`output`](@ref) but adding the result.
"""
lx_show(lxc::LxCom, lxd) = lx_output(lxc, lxd; res=true)

"""
    \\textinput{...}

Resolve a `\\textinput{rpath}` command: insert **and reprocess** some text.
If you just want to include text as a plaintext block use
`\\input{plaintext}{rpath}` instead.
"""
function lx_textinput(lxc::LxCom, lxd::Vector{LxDef})
    rpath = stent(lxc.braces[1])
    input = ""
    try
        fp,   = resolve_rpath(rpath)
        input = read(fp, String)
    catch e
        return html_err(e.m)
    end
    return reprocess(input, lxd)
end

"""
    \\figalt{...}{...}

Resolve a `\\figalt{alt}{rpath}` (find a fig and include it with alt).
"""
function lx_figalt(lxc::LxCom, _)
    rpath = stent(lxc.braces[2])
    alt   = stent(lxc.braces[1])
    path  = parse_rpath(rpath; canonical=false, code=true)
    fdir, fext = splitext(path)

    # there are several cases
    # A. a path with no extension --> guess extension
    # B. a path with extension --> use that
    # then in both cases there can be a relative path set but the user may mean
    # that it's in the subfolder /output/ (if generated by code) so should look
    # both in the relpath and if not found and if /output/ not already last dir
    candext = ifelse(isempty(fext),
                     (".png", ".jpeg", ".jpg", ".svg", ".gif", ".json"), (fext,))
    for ext ∈ candext
        candpath = fdir * ext
        syspath  = joinpath(PATHS[:site], split(candpath, '/')...)
        isfile(syspath) && return ext == ".json" ? html_plotly(candpath) : html_img(candpath, alt)
    end
    # now try in the output dir just in case (provided we weren't already
    # looking there)
    p1, p2 = splitdir(fdir)
    if splitdir(p1)[2] != "output"
        for ext ∈ candext
            candpath = joinpath(p1, "output", p2 * ext)
            syspath  = joinpath(PATHS[:site], split(candpath, '/')...)
            isfile(syspath) && return ext == ".json" ? html_plotly(candpath) : html_img(candpath, alt)
        end
    end
    return html_err("Image matching '$path' not found.")
end

"""
    \\tableinput{...}{...}

Resolve a `\\tableinput{header}{rpath}` (find a table+header and include it).
"""
function lx_tableinput(lxc::LxCom, _)
    rpath  = stent(lxc.braces[2])
    header = stent(lxc.braces[1])
    path   = parse_rpath(rpath; canonical=false)
    fdir, fext = splitext(path)
    # copy-paste from resolve_lx_figalt()
    # A. a path with extension --> use that
    # there can be a relative path set but the user may mean
    # that it's in the subfolder /output/ (if generated by code) so should look
    # both in the relpath and if not found and if /output/ not already the last subdir
    syspath = joinpath(PATHS[:site], split(path, '/')...)
    isfile(syspath) && return csv2html(syspath, header)
    # now try in the output dir just in case (provided we weren't already
    # looking there)
    p1, p2 = splitdir(fdir)
    if splitdir(p1)[2] != "output"
        candpath = joinpath(p1, "output", p2 * fext)
        syspath = joinpath(PATHS[:site], split(candpath, '/')...)
        isfile(syspath) && return csv2html(syspath, header)
    end
    return html_err("Table matching '$path' not found.")
end

"""
    \\literate{...}

Resolve a `\\literate{rpath}` (find a literate script and insert it).
"""
function lx_literate(lxc::LxCom, lxd::Vector{LxDef})
    rpath = stent(lxc.braces[1])
    opath, haschanged = literate_to_franklin(rpath)
    # check file is there
    if isempty(opath)
        return html_err("Literate file matching '$rpath' not found.")
    end
    if haschanged && FD_ENV[:FULL_PASS]
        set_var!(LOCAL_VARS, "reeval", true)
    end
    # then reprocess
    return reprocess(read(opath, String), lxd, nostripp=true)
end

"""
    \\nonumber{...}

Suppress the displaying of a number.
"""
function lx_nonumber(c, lxd)
    PAGE_EQREFS[PAGE_EQREFS_COUNTER] -= 1
    return reprocess("@@nonumber " * content(c.braces[1]) * "@@", lxd)
end