text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module MadNLPPardisoMKL

import ..MadNLP:
    @kwdef, Logger, @debug, @warn, @error,
    SubVector, StrideOneVector, SparseMatrixCSC, libblas,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    AbstractOptions, AbstractLinearSolver, set_options!,
    introduce, factorize!, solve!, improve!, is_inertia, inertia, blas_num_threads

const INPUT_MATRIX_TYPE = :csc

@enum(MatchingStrategy::Int,COMPLETE=1,COMPLETE2x2=2,CONSTRAINTS=3)

@kwdef mutable struct Options <: AbstractOptions
    pardisomkl_num_threads::Int = 1
    pardiso_matching_strategy::MatchingStrategy = COMPLETE2x2
    pardisomkl_max_iterative_refinement_steps::Int = 1
    pardisomkl_msglvl::Int = 0
    pardisomkl_order::Int = 2
end

mutable struct Solver <: AbstractLinearSolver
    pt::Vector{Ptr{Cvoid}}
    iparm::Vector{Int32}
    perm::Vector{Int32}
    msglvl::Ref{Int32}
    err::Ref{Int32}
    csc::SparseMatrixCSC{Float64,Int32}
    w::Vector{Float64}
    opt::Options
    logger::Logger
end

pardisomkl_pardisoinit(pt,mtype::Ref{Cint},iparm::Vector{Cint}) =
    ccall(
        (:pardisoinit,libblas),
        Cvoid,
        (Ptr{Cvoid},Ptr{Cint},Ptr{Cint}),
        pt,mtype,iparm)
pardisomkl_pardiso(pt,maxfct::Ref{Cint},mnum::Ref{Cint},mtype::Ref{Cint},
                   phase::Ref{Cint},n::Ref{Cint},a::Vector{Cdouble},ia::Vector{Cint},ja::Vector{Cint},
                   perm::Vector{Cint},nrhs::Ref{Cint},iparm::Vector{Cint},msglvl::Ref{Cint},
                   b::StrideOneVector{Cdouble},x::StrideOneVector{Cdouble},err::Ref{Cint}) =
                       ccall(
                           (:pardiso,libblas),
                           Cvoid,
                           (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
                            Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
                            Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},Ptr{Cint}),
                           pt,maxfct,mnum,mtype,phase,n,a,ia,ja,perm,nrhs,iparm,msglvl,b,x,err)

function pardisomkl_set_num_threads!(n)
    ccall((:mkl_set_dynamic, libblas),
         Cvoid,
         (Ptr{Cint},),
          Ref{Cint}(0))
    ccall((:mkl_set_num_threads, libblas),
          Cvoid,
          (Ptr{Cint},),
          Ref{Cint}(n))
end

function Solver(csc::SparseMatrixCSC;
                opt=Options(),logger=Logger(),
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                kwargs...)
    !isempty(kwargs) && (for (key,val) in kwargs; option_dict[key]=val; end)
    set_options!(opt,option_dict)

    w   = Vector{Float64}(undef,csc.n)

    pt = Vector{Ptr{Cvoid}}(undef,64)
    iparm = Vector{Int32}(undef,64)

    perm = Vector{Int32}(undef,csc.n)
    msglvl = Ref{Int32}(0)
    err = Ref{Int32}(0)

    pt.=0
    iparm[1]=0

    pardisomkl_pardisoinit(pt,Ref{Int32}(-2),iparm)

    iparm[1]=1
    iparm[2]=opt.pardisomkl_order # METIS
    iparm[3]=0
    iparm[6]=1
    iparm[8]=opt.pardisomkl_max_iterative_refinement_steps
    iparm[10]=12
    iparm[11]=2
    iparm[13]=1 # matching strateg
    iparm[21]=3 # bunch-kaufman pivotin
    iparm[24]=1 # parallel factorization
    iparm[25]=0 # parallel solv

    pardisomkl_set_num_threads!(opt.pardisomkl_num_threads)
    pardisomkl_pardiso(pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(11),
                     Ref{Int32}(csc.n),csc.nzval,csc.colptr,csc.rowval,perm,
                       Ref{Int32}(1),iparm,msglvl,Float64[],Float64[],err)
    pardisomkl_set_num_threads!(blas_num_threads[])

    err.x < 0  && throw(SymbolicException())

    M = Solver(pt,iparm,perm,msglvl,err,csc,w,opt,logger)

    finalizer(finalize,M)

    return M
end
function factorize!(M::Solver)
    pardisomkl_set_num_threads!(M.opt.pardisomkl_num_threads)
    pardisomkl_pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(22),
                     Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
                       Ref{Int32}(1),M.iparm,M.msglvl,Float64[],Float64[],M.err)
    M.err.x < 0  && throw(FactorizationException())
    pardisomkl_set_num_threads!(blas_num_threads[])

    M.iparm[14] == 0 && return M
    pardisomkl_set_num_threads!(M.opt.pardisomkl_num_threads)
    pardisomkl_pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(12),
                     Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
                       Ref{Int32}(1),M.iparm,M.msglvl,Float64[],Float64[],M.err)
    pardisomkl_set_num_threads!(blas_num_threads[])
    M.err.x < 0  && throw(FactorizationException())
    return M
end
function solve!(M::Solver,rhs::StrideOneVector{Float64})
    pardisomkl_set_num_threads!(M.opt.pardisomkl_num_threads)
    pardisomkl_pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(33),
                     Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
                       Ref{Int32}(1),M.iparm,M.msglvl,rhs,M.w,M.err)
    pardisomkl_set_num_threads!(blas_num_threads[])
    M.err.x < 0  && throw(SolveException())
    return rhs
end

function finalize(M::Solver)
    pardisomkl_pardiso(M.pt,Ref{Int32}(1),Ref{Int32}(1),Ref{Int32}(-2),Ref{Int32}(-1),
                     Ref{Int32}(M.csc.n),M.csc.nzval,M.csc.colptr,M.csc.rowval,M.perm,
                     Ref{Int32}(1),M.iparm,M.msglvl,Float64[],M.w,M.err)
end
is_inertia(::Solver)=true
function inertia(M::Solver)
    pos = M.iparm[22]
    neg = M.iparm[23]
    return (pos,M.csc.n-pos-neg,neg)
end

function improve!(M::Solver)
    @debug(M.logger,"improve quality failed.")
    return false
end
introduce(::Solver)="pardiso-mkl"

end # module


text:
TYPE_ID=1;
FAKE_ID=2;
DMO_ID=3;
VAR_ID=4;
FIX_ID=5;
DMO_ID=5;
_TYPE = 1;
_ID = 2;
_P = 3;
_PMAX = 4;
_μ = 6;
_σ = 5;

function read_id_fuel(dir_path::String)
    # id_fuel = read_id_fuel_txt(dir_path);
    id_fuel = read_id_fuel_csv(dir_path);
    return id_fuel;
end
function read_id_fuel_txt(dir_path::String)
    id_fuel = Dict{String, String}();
    file_path = joinpath(dir_path, "id_fuel.txt")
    open(file_path) do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = AmplTxt.split_with_space(ln);
                id_fuel[buffer[1]] = buffer[2];
            end
        end
    end
    return id_fuel;
end
function read_id_fuel_csv(dir_path::String)
    id_fuel = Dict{String, String}();
    file_path = joinpath(dir_path, "generators.csv")
    open(file_path) do file
        i = 0;
        for ln in eachline(file)
            i+=1;
            if i>1
                buffer = split(ln, ";");
                id_fuel[buffer[1]] = buffer[2];
            end
        end
    end
    return id_fuel;
end
function read_uncertainties(dir_path::String)
    uncertainties = Vector{String}[];
    pmax_uncertain = Dict{String, Float64}();
    file_path = joinpath(dir_path, "capa_wind.txt");
    open(file_path) do file
        for ln in eachline(file)
            # don't read commentted line
            if ln[1] != '#'
                buffer = AmplTxt.split_with_space(ln);
                push!(uncertainties, buffer)
                pmax_uncertain[buffer[_ID]] = parse(Float64, buffer[_PMAX]);
            end
        end
    end
    return uncertainties, pmax_uncertain;
end

function read_bus_load(ampltxt)
    buses = amplTxt["buses"];
    num_to_name= Dict{Int64, String}()
    bus_load = Dict{Int64, Float64}()
    for bus in buses.data
        num = parse(Int, bus[2]);
        name = bus[11];
        num_to_name[num] = name;
        bus_load[num]  = 0;
    end
    loads = amplTxt["loads"];
    for load in loads.data
        num = parse(Int, load[2]);
        bus = parse(Int, load[3]);
        p = parse(Float64, load[11]);
        bus_load[bus] = get(bus_load, bus, 0) + p;
    end
    return num_to_name, bus_load;

end

function write_pscopf_previsions(output_file_path, uncertainties, SCENARIO, HORIZON, TIME_STEP)
    open(output_file_path, "w") do file
        write(file, @sprintf("#%-9s%25s%25s%10s\n", "id", "h_15m", "ech", "v"))

        generators = amplTxt["generators"];
        for genData in generators.data
            gen = parse(Int, genData[2]);
            p = parse(Float64, genData[16]);
            name = genData[20];
            for ech in HORIZON, ts in TIME_STEP
                real_h = h+ts
                write(file, @sprintf("%-10s%25s%25s%10.3f\n", name,real_h, ech,  p))
            end
        end
    end
end

function write_pscopf_uncertainties(output_file_path, uncertainties, SCENARIO, HORIZON, TIME_STEP)
    open(output_file_path, "w") do file
        write(file, @sprintf("#%-9s%25s%25s%10s%10s\n", "id", "h_15m", "ech","scenario", "v"))
        is_uncertain = Dict{String, Bool}();
        for uncertainty in uncertainties
            id = uncertainty[_ID];
            is_uncertain[id] = true;
            for ech in HORIZON, ts in TIME_STEP
                p = parse(Float64, uncertainty[_P]);
                pMax = parse(Float64, uncertainty[_PMAX]);
                real_h = h+ts
                factor=1
                if ech == h
                    factor = 0.5
                elseif ech==h_2
                    factor = 2
                end

                μ = factor*parse(Float64, uncertainty[_μ]);
                σ = factor*parse(Float64, uncertainty[_σ]);

                normal = rand(Normal(μ, σ), N_SCENARIO);
                i = 0;
                for scenario in SCENARIO
                    i += 1
                    v =  p + normal[i] * pMax
                    v = min(pMax, v)
                    v = max(0, v)
                    # println(@sprintf("%10s%10.3f%10.3f", id, p, v))
                    write(file, @sprintf("%-10s%25s%25s%10s%10.3f\n", id, real_h,ech, scenario,  v))
                end
            end
        end
        generators = amplTxt["generators"];
        for genData in generators.data
            gen = parse(Int, genData[2]);
            p = parse(Float64, genData[16]);
            name = genData[20];
            if !haskey(is_uncertain, name)
                for ech in HORIZON, scenario in SCENARIO, ts in TIME_STEP
                    real_h = h+ts
                    write(file, @sprintf("%-10s%25s%25s%10s%10.3f\n", name,real_h, ech, scenario,  p))
                end
            end
        end

        for kvp in bus_load
            bus = kvp[1]
            name = num_to_name[bus]
            v = kvp[2]
            for ech in HORIZON, scenario in SCENARIO, ts in TIME_STEP
                real_h = h+ts
                write(file,  @sprintf("%-10s%25s%25s%10s%10.3f\n",name,real_h, ech, scenario, v))
            end
        end

    end
end

function write_pscopf_units(output_file_path, amplTxt, id_fuel, flexibility)
    open(output_file_path, "w") do file
        write(file, @sprintf("#%-9s %13s %13s %13s %13s %13s %13s\n", "name", "p", "minP", "maxP", "start", "prop", "dmo"))
        generators = amplTxt["generators"];
        for genData in generators.data
            gen = parse(Int, genData[2]);
            bus = parse(Int, genData[3]);
            conbus = parse(Int, genData[4]);
            minP = parse(Float64, genData[6])+10;
            maxP = parse(Float64, genData[7]);
            p = parse(Float64, genData[16]);
            name = genData[20];

            fuel = id_fuel[name]

            carac = flexibility[fuel]

            start_cost = carac[FIX_ID];
            prop_cost = carac[VAR_ID];
            dmo = carac[DMO_ID];
            write(file, @sprintf("%-10s %.8E %.8E %.8E %.8E %.8E %.8E\n", name, p, minP, maxP, start_cost, prop_cost, dmo))
        end
    end

end

function write_pscopf_gen_type_bus(output_file_path, amplTxt, id_fuel, num_to_name)
    open(output_file_path, "w") do file
        generators = amplTxt["generators"];
        for genData in generators.data
            gen = parse(Int, genData[2]);
            bus = parse(Int, genData[3]);
            conbus = parse(Int, genData[4]);
            name = genData[20];

            fuel = id_fuel[name]
            if fuel in ["onshore", "pv"]
                kind  ="Limitable"
            else
                kind = "Imposable"
            end
            write(file, @sprintf("%s %s %s\n",name, kind,  num_to_name[conbus]))
        end
    end

end

function read_flexibility(dir_path)
    # flexibility =read_flexibility_hard(dir_path);
    flexibility =read_flexibility_csv(dir_path);
    return flexibility;
end

function read_flexibility_hard(dir_path)
    flexibility = Dict([
        "coal"=>["coal", 600, 360, 50, 180000],
        "nuclear"=>["nuclear", 1000, 1440, 10, 1000000],
        "gas_conventional"=>["gas", 200, 300, 40, 40000],
        "gas_ccg"=>["gas", 500, 180, 40, 40000],
        "gas_ct"=>["gas", 200, 30, 120, 12000],
        "onshore"=>["wind", 0, 5, 50, 0],
        "pv"=>["solar", 0, 5, 50, 0]
    ]);
end
function read_flexibility_csv(dir_path)
    flexibility = Dict{String, Vector{Any}}();
    file_path = joinpath(dir_path, "flexibility.csv")
    open(file_path) do file
        i = 0;
        for ln in eachline(file)
            i+=1;
            if i>2
                buffer = split(ln, ";");
                name = buffer[1];
                type = buffer[2];
                p = parse(Float64,buffer[3]);
                dmo =parse(Float64,buffer[4]);
                variable =parse(Float64,buffer[5]);
                fixed  =parse(Float64,buffer[6]);
                push!(flexibility, type=>[name, p, dmo, variable, fixed]);
            end
        end
    end
    return flexibility;
end

ext: jl

max_stars_repo_name: azized/pscopf

max_issues_repo_name: azized/pscopf

max_forks_repo_name: azized/pscopf


text:
struct MyStruct
    a
    c
end

"""
overide the print version
"""
function show(io::IO, m::MyStruct)
    print("2 arg version of MyStuct - we be used with print")
end

"""
overide the representation method
"""
function show(io::IO, ::MIME"text/plain",  m::MyStruct)
    print("3 arg version of MyStuct - excuted when called bare")
end

module Testy
using Random

export a

function a()
    println("this is a")
end

# Top-Level code will be excuted when including / using
println ("this is a top level code")

end # module Testy

# We can't create and use a new version of the Module because it will create a conflict which version of
# the new method to overide
#

module Testy
using Random

export a

function a()
    println("this is a")
    println("changed a")
end

# Top-Level code will be excuted when including / using
println ("this is a top level code")

end # module Testy


text:
# Bogumił Kamiński, 2022

# Codes for chapter 13

# Codes for section 13.1

using CSV
using CategoricalArrays
using DataFrames
using DataFramesMeta
using Dates
using Distributions
import Downloads
using FreqTables
using GLM
using Plots
using Random
using ROCAnalysis
using SHA
using Statistics
import ZipFile

url_zip = "https://stacks.stanford.edu/file/druid:yg821jf8611/" *
          "yg821jf8611_ky_owensboro_2020_04_01.csv.zip";
local_zip = "owensboro.zip";
isfile(local_zip) || Downloads.download(url_zip, local_zip)
isfile(local_zip)
open(sha256, local_zip) == [0x14, 0x3b, 0x7d, 0x74,
                            0xbc, 0x15, 0x74, 0xc5,
                            0xf8, 0x42, 0xe0, 0x3f,
                            0x8f, 0x08, 0x88, 0xd5,
                            0xe2, 0xa8, 0x13, 0x24,
                            0xfd, 0x4e, 0xab, 0xde,
                            0x02, 0x89, 0xdd, 0x74,
                            0x3c, 0xb3, 0x5d, 0x56]
archive = ZipFile.Reader(local_zip)
owensboro = @chain archive begin
    only(_.files)
    read
    CSV.read(DataFrame; missingstring="NA")
end;
close(archive)

summary(owensboro)
describe(owensboro, :nunique, :nmissing, :eltype)

# Code for listing 13.1

select!(owensboro, :date, :type, :arrest_made, :violation);
summary(owensboro)
describe(owensboro, :nunique, :nmissing, :eltype)

# Code for listing 13.2

df = DataFrame(id=[1, 2, 1, 2], v=1:4)
combine(df, :v => sum => :sum)
transform(df, :v => sum => :sum)
select(df, :v => sum => :sum)
gdf = groupby(df, :id)
combine(gdf, :v => sum => :sum)
transform(gdf, :v => sum => :sum)
select(gdf, :v => sum => :sum)

# Code for listing 13.3

select(df,
    :v => identity => :v1,
    :v => identity,
    :v => :v2,
    :v)

# Codes for section 13.2

owensboro.violation

# Code for listing 13.4

violation_list = [strip.(split(x, ";"))
                  for x in owensboro.violation]
violation_flat = reduce(vcat, violation_list)
violation_flat_clean = [contains(x, "SPEEDING") ?
                        "SPEEDING" : x for x in violation_flat]
sort(freqtable(violation_flat_clean), rev=true)

# Code for listing 13.5

agg_violation = @chain owensboro begin
    select(:violation =>
           ByRow(x -> strip.(split(x, ";"))) =>
           :v)
    flatten(:v)
    select(:v =>
           ByRow(x -> contains(x, "SPEEDING") ? "SPEEDING" : x) =>
           :v)
    groupby(:v)
    combine(nrow => :count)
    sort(:count, rev=true)
end

# Code explaining ByRow

sqrt(4)
sqrt([4, 9, 16])
ByRow(sqrt)([4, 9, 16])
f = ByRow(sqrt)
f([4, 9, 16])

# Code explaining flatten

df = DataFrame(id=1:2, v=[[11, 12], [13, 14, 15]])
flatten(df, :v)

# Code explaining nrow

@chain DataFrame(id=[1, 1, 2, 2, 2]) begin
    groupby(:id)
    combine(nrow, nrow => :rows)
end

# Code explaining sort

df = DataFrame(a=[2, 1, 2, 1, 2], b=5:-1:1)
sort(df, :b)
sort(df, [:a, :b])

# Code showing an example transformation

df = DataFrame(x=[4, 9, 16])
transform(df, :x => ByRow(sqrt))

# Code showing usage of DataFramesMeta.jl

@chain owensboro begin
    @rselect(:v=strip.(split(:violation, ";")))
    flatten(:v)
    @rselect(:v=contains(:v, "SPEEDING") ?
                         "SPEEDING" : :v)
    groupby(:v)
    combine(nrow => :count)
    sort(:count, rev=true)
end

# Code showing comparison of DataFramesMeta.jl vs DataFrames.jl

df = DataFrame(x=[4, 9, 16])
@select(df, :s = sqrt.(:x))
@rselect(df, :s = sqrt(:x))
select(df, :x => ByRow(sqrt) => :s)

# Codes for section 13.3

# Code for listing 13.6

owensboro2 = select(owensboro,
    :arrest_made => :arrest,
    :date => ByRow(dayofweek) => :day,
    :type,
    [:violation =>
     ByRow(x -> contains(x, agg_violation.v[i])) =>
     "v$i" for i in 1:4])

# Code explaining programmatic generation of transformations

[:violation =>
 ByRow(x -> contains(x, agg_violation.v[i])) =>
 "v$i" for i in 1:4]

combine(owensboro, [:date :arrest_made] .=> [minimum, maximum])
[:date :arrest_made] .=> [minimum, maximum]

# Code for listing 13.7

weekdays = DataFrame(day=1:7,
                     dayname=categorical(dayname.(1:7); ordered=true))
isordered(weekdays.dayname)
levels(weekdays.dayname)
levels!(weekdays.dayname, weekdays.dayname)

# Code showing join operation

leftjoin!(owensboro2, weekdays; on=:day)

# Code for listing 13.8

@chain owensboro2 begin
    groupby([:day, :dayname]; sort=true)
    combine(nrow)
end

# Code showing creation of frequency table from a data frame

freqtable(owensboro2, :dayname, :day)

# Code for listing 13.9

@chain owensboro2 begin
    groupby([:day, :dayname]; sort=true)
    combine(nrow)
    unstack(:dayname, :day, :nrow; fill=0)
end

# Code for dropping rows with missing data

dropmissing!(owensboro2)

# Code for dropping unwanted column

select!(owensboro2, Not(:day))

# Codes for section 13.4

# Code for listing 13.10

Random.seed!(1234);
owensboro2.train = rand(Bernoulli(0.7), nrow(owensboro2));
mean(owensboro2.train)

# Code for listing 13.11

train = subset(owensboro2, :train)
test = subset(owensboro2, :train => ByRow(!))

# Code showing subsetting with DataFramesMeta.jl

@rsubset(owensboro2, !(:train))

# Code building a predictive model

model = glm(@formula(arrest~dayname+type+v1+v2+v3+v4),
            train, Binomial(), LogitLink())

# Code for making predictions

train.predict = predict(model)
test.predict = predict(model, test)

# Code for producing histograms showing predictions

test_groups = groupby(test, :arrest);
histogram(test_groups[(false,)].predict;
          bins=10, normalize=:probability,
          fillstyle= :/, label="false")
histogram!(test_groups[(true,)].predict;
           bins=10, normalize=:probability,
           fillalpha=0.5, label="true")

# Code for computing of confusion matrix

@chain test begin
    @rselect(:predicted=:predict > 0.15, :observed=:arrest)
    proptable(:predicted, :observed; margins=2)
end

# Code for listing 13.12

test_roc =  roc(test; score=:predict, target=:arrest)
plot(test_roc.pfa, test_roc.pmiss;
     color="black", lw=3,
     label="test (AUC=$(round(100*auc(test_roc), digits=2))%)",
     xlabel="pfa", ylabel="pmiss")
train_roc =  roc(train, score=:predict, target=:arrest)
plot!(train_roc.pfa, train_roc.pmiss;
      color="gold", lw=3,
      label="train (AUC=$(round(100*auc(train_roc), digits=2))%)")

# Codes for section 13.5

DataFrame(:a=>1, :a=>2)

DataFrame(:a=>1, :a=>2; makeunique=true)

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
@testset "LTI system kernels tests" begin
    N = 10
    types = [Int32, Int64, Float32, Float64, ComplexF32, ComplexF64]
    Δt = 0.1
    atol = Δt

    @testset "Diff tests" begin

        @testset "identity operator Diff(0) tests" begin
            id = Romeo.LTI.Diff(0)
            k = Romeo.LTI.kernel(id, N, Δt)
            @test k[1] == 1 || "The first element of the identity kernel should be equal to 1"
            for i = 2:N
                @test k[i] == 0 || "All elements with index ≥2 of the identity kernel should be zero."
            end # for

            ys = Romeo.LTI.step_resp(id, N, Δt)
            for i = 1:N
                @test ys[i] == 1 || "All elements of the step response of the identity operator should be equal to 1."
            end # for
        end # testset

        @testset "first order integral operator tests" begin
            int = Romeo.LTI.Diff(-1)
            k = Romeo.LTI.kernel(int, N, Δt)
            for i = 1:N
                @test k[i] == Δt || "All elements of the first order integral kernel should be equal to Δt."
            end # for

            yi = Romeo.LTI.impulse_resp(int, N, Δt)
            for i = 1:N
                @test yi[i] == 1 || "All elements of the impulse response of the first order integrator should be equal to 1."
            end # for

            ys = Romeo.LTI.step_resp(int, N, Δt)
            for i = 1:N
                expected = (i - 1) * Δt
                @test isapprox(ys[i], expected; atol=atol) || "Step response of the first order integrator at index $i : $(ys[i]) ≠ $expected"
            end # for
        end # testset

        @testset "second order integral operator tests" begin

            int2 = Romeo.LTI.Diff(-2)
            k = Romeo.LTI.kernel(int2, N, Δt)
            for i = 1:N
                expected = 1 // 2 * Δt^2
                @test isapprox(k[i], expected; atol=atol) || "Element of the second order integral kernel at index $i : $(k[i]) ≠ $expected."
            end # for

        end # testset

        @testset "differentiator tests" begin
            Δt = 0.1
            s = Romeo.LTI.Diff{Float64}(1)

            k = Romeo.LTI.kernel(s, N, Δt)
            @test k[1] == 1 / Δt || "The first element of the differentiation kernel should be equal to 1/Δt"
            @test k[2] == -1 / Δt || "The second element of the differentiation kernel should be -1/Δt"
            for i = 3:N
                @test k[i] == 0 || "All elements with index ≥3 of the differentiation kernel should be zero."
            end # for

            ys = Romeo.LTI.step_resp(s, N, Δt)
            @test ys[1] == 1 / Δt || "The first element of the step response of a differentiator should be 1/Δt."
            for i = 2:N
                @test ys[i] == 0 || "All elements of the step response of a differentiator should be zero, except the first one."
            end # for
        end # testset

    end

end # testset


text:
import Base: @get!

include("read.jl")
include("dump.jl")
include("sugar.jl")

# Display

syntax(v::Vertex) = syntax(dl(v))

function Base.show(io::IO, v::Vertex)
  print(io, typeof(v))
  print(io, "(")
  s = MacroTools.alias_gensyms(syntax(v))
  if length(s.args) == 1
    print(io, sprint(print, s.args[1]))
  else
    foreach(x -> (println(io); print(io, sprint(print, x))), s.args)
  end
  print(io, ")")
end

import Juno: Row, Tree

code(x) = Juno.Model(Dict(:type=>"code",:text=>x))

@render Juno.Inline v::Vertex begin
  s = MacroTools.alias_gensyms(syntax(v))
  Tree(typeof(v), map(s -> code(string(s)), s.args))
end

# Function / expression macros

export @flow, @iflow, @vtx, @ivtx

function inputsm(args)
  bindings = d()
  for arg in args
    isa(arg, Symbol) || error("invalid argument $arg")
    bindings[arg] = constant(arg)
  end
  return bindings
end

type SyntaxGraph
  args::Vector{Symbol}
  output::DVertex{Any}
end

function flow_func(ex)
  @capture(shortdef(ex), name_(args__) = exs__)
  bs = inputsm(args)
  output = graphm(bs, exs)
  :($(esc(name)) = $(SyntaxGraph(args, output)))
end

function flowm(ex, f = dl)
  isdef(ex) && return flow_func(ex)
  g = graphm(block(ex))
  g = mapconst(x -> isexpr(x, :$) ? esc(x.args[1]) : Expr(:quote, x), g)
  constructor(f(g))
end

macro flow(ex)
  flowm(ex)
end

macro iflow(ex)
  flowm(ex, il)
end

function vtxm(ex, f = dl)
  exs = graphm(block(ex))
  @>> exs graphm mapconst(esc) f constructor
end

macro vtx(ex)
  vtxm(ex)
end

macro ivtx(ex)
  vtxm(ex, il)
end


text:
module find

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using GeoData, DimensionalData, Test, Statistics, Dates
using GeoData.LookupArrays, GeoData.Dimensions

data1 = cumsum(cumsum(ones(10, 11); dims=1); dims=2)
data2 = 2cumsum(cumsum(ones(10, 11, 1); dims=1); dims=2)
dims1 = X(10.0:10.0:100.0), Y(-50.0:10.0:50.0)
dims2 = (dims1..., Ti([DateTime(2019)]))
refdimz = ()
nme = :test
mval = -9999.0
meta = NoMetadata()

# Formatting only occurs in shorthand constructors
ga1 = GeoArray(data1, dims1; refdims=refdimz, name=nme, metadata=meta, missingval=mval)
ga2 = GeoArray(data2, dims2)

@testset "stack constructors" begin
    st1 = GeoStack((ga1, ga2); name=(:ga1, :ga2))
    st2 = GeoStack((ga1, ga2); keys=(:ga1, :ga2))
    st3 = GeoStack((data1, data2), dims2; keys=(:ga1, :ga2))
    st4 = GeoStack(st3)
    st5 = GeoStack(st3, dims2)
    @test st1 == st2 == st3 == st4 == st5

    # The dimension differences are lost because the table
    # is tidy - every column is the same length
    table_st = GeoStack(DimTable(st3), dims2)
    @test dims(table_st[:ga1]) isa Tuple{<:X,<:Y,<:Ti}
end

st = GeoStack((ga1, ga2); name=(:ga1, :ga2))

@testset "stack layers" begin
    @test length(st) == 2
    @test first(st) == ga1
    @test last(st) == ga2
    @test DimensionalData.layers(st) isa NamedTuple
    @test st[:ga1] == ga1
    @test st[:ga2] == ga2
    @test parent(st[:ga1]) == data1
    @test parent(st[:ga1]) isa Array{Float64,2}
    @test keys(st) == (:ga1, :ga2)
    @test haskey(st, :ga1)
    @test names(st) == (:ga1, :ga2)
    @test collect(values(st)) == [ga1, ga2]
    # We can iterate/splat stacks as GeArrays
    @test st == GeoStack(st...)
end

@testset "st fields " begin
    @test DimensionalData.layerdims(st, :ga1) == DimensionalData.format(dims1, data1)
    @test metadata(st) == NoMetadata()
    @test metadata(st, :ga1) == NoMetadata()
end

@testset "indexing" begin
    # Indexing the st is the same as indexing its child array
    a = st[:ga1][X(2:4), Y(5:6)]
    @test a == st[:ga1, X(2:4), Y(5:6)]

    @inferred st[:ga1][X=2:4, Y=5:6]
    # FIXME: This isn't inferred, the constants don't propagate like they
    # do in the above call. Probably due to the anonymous wrapper function.
    @test_broken @inferred st[:ga1, X(2:4), Y(5:6)] isa GeoArray

    # Getindex for a whole st of new GeoArrays
    a = st[X=2:4, Y=5:6]
    @test a isa GeoStack
    @test a[:ga1] isa GeoArray
    @test parent(a[:ga1]) isa Array
    @test a[:ga1] == data1[2:4, 5:6]
    @test a[:ga2] == data2[2:4, 5:6, 1:1]

    @testset "select new arrays for the whole st" begin
        s = st[Y=Between(-10, 10.0), Ti=At(DateTime(2019))]
        @test s isa GeoStack
        @test s[:ga1] isa GeoArray
        @test parent(s[:ga1]) isa Array
        @test s[:ga1] == data1[:, 5:7]
        @test s[:ga2] == data2[:, 5:7, 1]
        @test dims(s[:ga2]) == (X(Sampled(10.0:10.0:100.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())),
                                Y(Sampled(-10.0:10.0:10.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())))
        @test refdims(s[:ga2]) == (Ti(Sampled([DateTime(2019)], ForwardOrdered(), Irregular(nothing, nothing), Points(), NoMetadata())),)
        @test ismissing(missingval(s, :ga2)) && ismissing(missingval(s[:ga2]))
    end

    @testset "select views of arrays for the whole st" begin
        sv = view(st, Y=Between(-4.0, 27.0), Ti=At(DateTime(2019)))
        @test sv isa GeoStack
        @test sv[:ga1] isa GeoArray
        @test parent(sv[:ga1]) isa SubArray
        @test sv[:ga1] == data1[:, 6:8]
        @test sv[:ga2] == data2[:, 6:8, 1]
        @test dims(sv[:ga2]) == (X(Sampled(10.0:10:100.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())),
                                 Y(Sampled(0.0:10:20.0, ForwardOrdered(), Regular(10.0), Points(), NoMetadata())))
        @test refdims(sv[:ga2])[1] == Ti(Sampled(view([DateTime(2019)], 1:1), ForwardOrdered(), Irregular((nothing, nothing)), Points(), NoMetadata()))
        # Stack of view-based GeoArrays
        v = view(st, X(2:4), Y(5:6))
        # @inferred view(st, X(2:4), Y(5:6))
        @test v isa GeoStack
        @test v[:ga1] isa GeoArray
        @test parent(v[:ga1]) isa SubArray
        @test v[:ga1] == view(data1, 2:4, 5:6)
        @test v[:ga2] == view(data2, 2:4, 5:6, 1:1)
    end

end

@testset "subset st with specific key(s)" begin
    s1 = GeoStack(st; keys=(:ga2,))
    @test keys(s1) == (:ga2,)
    @test length(values(s1)) == 1
    s2 = GeoStack(st; keys=(:ga1, :ga2))
    @test keys(s2) == (:ga1, :ga2)
    @test length(values(s2)) == 2
end

@testset "concatenate stacks" begin
    dims1b = X(110:10:200), Y(-50:10:50)
    dims2b = (dims1b..., Ti([DateTime(2019)]))
    stack_a = GeoStack((ga1=ga1, ga2=ga2))
    stack_b = GeoStack((ga1=GeoArray(data1 .+ 10, dims1b), ga2=GeoArray(data2 .+ 20, dims2b)))
    catstack = cat(stack_a, stack_b; dims=X())
    @test size(first(catstack)) == (20, 11)
    @test val(dims(catstack, X)) ≈ 10.0:10.0:200.0
    #@test step(dims(first(catstack), X())) == 10.0
    @test DimensionalData.bounds(dims(first(catstack), X)) == (10.0, 200.0)
    @test catstack[:ga1][Y(1)] == 1.0:20.0
    @test catstack[:ga2][Y(1), Ti(1)] == 2.0:2.0:40.0
    catstack = cat(stack_a, stack_b; dims=(X(), Y()))
    @test size(first(catstack)) == (20, 22)
end

@testset "copy" begin
    cp = copy(st)
    @test all(st[:ga1] .=== cp[:ga1])
    @test st[:ga1] !== cp[:ga1]
end

@testset "show" begin
    sh = sprint(show, st)
    # Test but don't lock this down too much
    @test occursin("GeoStack", sh)
    @test occursin("Y", sh)
    @test occursin("X", sh)
end


text:
# indicator of a point

export IndPoint

"""
**Indicator of a singleton**

    IndPoint(p=0.0)

Returns the indicator of the set
```math
C = \\{p \\}.
```
Parameter `p` can be a scalar, in which case the unique element of `S` has uniform coefficients.
"""
struct IndPoint{T <: Union{Real, Complex, AbstractArray{<:RealOrComplex}}} <: ProximableFunction
  p::T
  function IndPoint{T}(p::T) where {T <: Union{Real, Complex, AbstractArray{<:RealOrComplex}}}
    new(p)
  end
end

is_separable(f::IndPoint) = true
is_convex(f::IndPoint) = true
is_singleton(f::IndPoint) = true
is_cone(f::IndPoint) = norm(f.p) == 0
is_affine(f::IndPoint) = true

IndPoint(p::T=0.0) where {T <: Union{Real, Complex, AbstractArray{<:RealOrComplex}}} = IndPoint{T}(p)

function (f::IndPoint{R})(x::AbstractArray{T}) where {R <: RealOrComplex, T <: RealOrComplex}
  for k in eachindex(x)
    if abs(x[k] - f.p) > 1e-14
      return +Inf
    end
  end
  return zero(T)
end

function (f::IndPoint{A})(x::AbstractArray{T}) where {A <: AbstractArray, T <: RealOrComplex}
  for k in eachindex(x)
    if abs(x[k] - f.p[k]) > 1e-14
      return +Inf
    end
  end
  return zero(T)
end

function prox!(y::AbstractArray{T}, f::IndPoint{R}, x::AbstractArray{T}, gamma::Real=1.0) where {R <: RealOrComplex, T <: RealOrComplex}
  for k in eachindex(x)
    y[k] = f.p
  end
  return zero(T)
end

function prox!(y::AbstractArray{T}, f::IndPoint{A}, x::AbstractArray{T}, gamma::Real=1.0) where {A <: AbstractArray, T <: RealOrComplex}
  for k in eachindex(x)
    y[k] = f.p[k]
  end
  return zero(T)
end

prox!(y::AbstractArray{T}, f::IndPoint, x::AbstractArray{T}, gamma::AbstractArray) where {T <: RealOrComplex} = prox!(y, f, x, 1.0)

fun_name(f::IndPoint) = "indicator of a point"
fun_dom(f::IndPoint) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndPoint) = "x ↦ 0 if x = p, +∞ otherwise"
fun_params(f::IndPoint) =
  string( "p = ", typeof(f.p) <: AbstractArray ? string(typeof(f.p), " of size ", size(f.p)) : f.p, ", ")

function prox_naive(f::IndPoint, x::AbstractArray{T}, gamma=1.0) where T <: RealOrComplex
    y = similar(x)
    y .= f.p
    return y, 0.0
end


text:
module ManifoldDiffEq

using ManifoldsBase
using Manifolds
using SciMLBase:
    SciMLBase, AbstractODEProblem, AbstractODEFunction, NullParameters, promote_tspan
import SciMLBase: build_solution
using OrdinaryDiffEq

using OrdinaryDiffEq:
    InterpolationData,
    OrdinaryDiffEqAlgorithm,
    OrdinaryDiffEqAdaptiveAlgorithm,
    OrdinaryDiffEqMutableCache,
    OrdinaryDiffEqConstantCache,
    trivial_limiter!,
    constvalue,
    @muladd,
    @unpack,
    @cache,
    @..

import OrdinaryDiffEq: alg_cache, alg_order, initialize!, perform_step!

using RecursiveArrayTools

include("utils.jl")
include("error_estimation.jl")
include("operators.jl")
include("problems.jl")
include("interpolation.jl")
include("frozen_solvers.jl")
include("lie_solvers.jl")

end # module


text:
using RayTracer
using Test

@testset "Tuples" begin
    include("tuples.jl")
end

@testset "Canvas" begin
    include("canvas.jl")
end

@testset "Transformations" begin
    include("transformations.jl")
end

@testset "Rays" begin
    include("rays.jl")
end

@testset "Sphere" begin
    include("sphere.jl")
end

@testset "Light" begin
    include("light.jl")
end

@testset "Materials" begin
    include("material.jl")
end


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

cxxinclude("TBufferJSON.h")

export rootjson

rootjson(obj::CppValue, compact::Integer = 3) =
	rootjson(pointer_to(obj), compact)

rootjson(obj::CppPtr, compact::Integer = 3) =
	string(@cxx TBufferJSON::ConvertToJSON(obj, Int32(compact)))


text:
quadratic_field(10)
k, a = cyclotomic_field(11)
automorphism_group(k)
kt, t = k["t"]
factor(t^2-a)
factor(t^5-a)
k, a = wildanger_field(3, 13)
h = hilbert_class_field(k)
K = number_field(h)
L = simple_extension(K)[1]
absolute_simple_field(L)
discriminant(maximal_order(K))
norm_equation(k, 27)


text:
macro Ax_max()
	code = quote
		# maximum longitudinal acceleration for given speed
		Ax_max = zeros(length(U),1)
		for i in eachindex(U)
			Ax_max[i,1] = AXC[1]*U[i]^3 + AXC[2]*U[i]^2 + AXC[3]*U[i] + AXC[4]
		end
		Ax_max
	end
	return esc(code)
end


text:
export @Eval

macro Eval(expr...)
	for ex in expr
	eval(ex)
	end
end


text:
#--------------------------------------------------------------------#
#     Input data used for heat and moisture transfer calculation     #
#--------------------------------------------------------------------#

#DODATI RUBNE UVJETE

#########
#Timestep if necessary#
##########

dt=60*60.        #Timestep in h
T=2919.      #Time domain

###########
#Mesh size#
###########

dl=0.005 #defined in m

########################
#Geometry and materials#
########################

#Elements -- matrix of vectors E_1=[material_key,x_0,y_0,x_1,y_1] to E_n

ne=27 #number of elements

E=zeros(Float64,ne,9) #Geometry and material matrix + boundary conditions

#Element 1 (Concrete)
E[1,1]=12
E[1,2]=0
E[1,3]=0
E[1,4]=0.76
E[1,5]=0.08
E[1,6]=1
E[1,7]=1.65
E[1,8]=2200
E[1,9]=1000

#Element 2 (XPS)
E[2,1]=14
E[2,2]=0
E[2,3]=0.08
E[2,4]=0.08
E[2,5]=1
E[2,6]=2
E[2,7]=0.035
E[2,8]=50
E[2,9]=1015

#Element 3 (Reinforced concrete)
E[3,1]=13
E[3,2]=0.08
E[3,3]=0.08
E[3,4]=0.68
E[3,5]=1
E[3,6]=3
E[3,7]=2.6
E[3,8]=2500
E[3,9]=1000

#Element 4 (Reinforced concrete)
E[4,1]=13
E[4,2]=0.16
E[4,3]=1
E[4,4]=0.68
E[4,5]=1.08
E[4,6]=4
E[4,7]=2.6
E[4,8]=2500
E[4,9]=1000

#Element 5 (XPS)
E[5,1]=14
E[5,2]=0
E[5,3]=1
E[5,4]=0.08
E[5,5]=1.08
E[5,6]=5
E[5,7]=0.035
E[5,8]=50
E[5,9]=1015

#Element 6 (XPS)
E[6,1]=14
E[6,2]=0.08
E[6,3]=1.01
E[6,4]=0.15
E[6,5]=1.08
E[6,6]=6
E[6,7]=0.035
E[6,8]=50
E[6,9]=1015

#Element 7 (PVC hydroinsulation)
E[7,1]=15
E[7,2]=0.08
E[7,3]=1
E[7,4]=0.16
E[7,5]=1.01
E[7,6]=7
E[7,7]=0.14
E[7,8]=1200
E[7,9]=1000

#Element 8 (PVC hydroinsulation)
E[8,1]=15
E[8,2]=0.15
E[8,3]=1.01
E[8,4]=0.16
E[8,5]=1.85
E[8,6]=8
E[8,7]=0.14
E[8,8]=1200
E[8,9]=1000

#Element 9 (XPS)
E[9,1]=14
E[9,2]=0.68
E[9,3]=0.08
E[9,4]=0.76
E[9,5]=1.08
E[9,6]=9
E[9,7]=0.035
E[9,8]=50
E[9,9]=1015

#Element 10 (Concrete)
E[10,1]=12
E[10,2]=0.76
E[10,3]=1
E[10,4]=1.4
E[10,5]=1.08
E[10,6]=10
E[10,7]=1.65
E[10,8]=2200
E[10,9]=1000

#Element 11 (XPS)
E[11,1]=14
E[11,2]=0.07
E[11,3]=1.08
E[11,4]=0.15
E[11,5]=1.85
E[11,6]=11
E[11,7]=0.035
E[11,8]=50
E[11,9]=1015

#Element 12 (Reinforced concrete)
E[12,1]=13
E[12,2]=0.16
E[12,3]=1.08
E[12,4]=0.36
E[12,5]=1.85
E[12,6]=12
E[12,7]=2.6
E[12,8]=2500
E[12,9]=1000

#Element 13 (PVC hydroinsulation)
E[13,1]=15
E[13,2]=0.36
E[13,3]=1.08
E[13,4]=0.37
E[13,5]=1.27
E[13,6]=13
E[13,7]=0.14
E[13,8]=1200
E[13,9]=1000

#Element 14 (PVC hydroinsulation)
E[14,1]=15
E[14,2]=0.37
E[14,3]=1.08
E[14,4]=1.4
E[14,5]=1.09
E[14,6]=14
E[14,7]=0.14
E[14,8]=1200
E[14,9]=1000

#Element 15 (Mineral wool)
E[15,1]=11
E[15,2]=0.37
E[15,3]=1.09
E[15,4]=1.40
E[15,5]=1.19
E[15,6]=15
E[15,7]=0.035
E[15,8]=150
E[15,9]=1030

#Element 16 (Cement screed)
E[16,1]=16
E[16,2]=0.37
E[16,3]=1.19
E[16,4]=1.4
E[16,5]=1.27
E[16,6]=16
E[16,7]=1.6
E[16,8]=2000
E[16,9]=1100

##############################
#Boundary condition materials#
##############################

#Element 17 (Constant temperature soil)
E[17,1]=100
E[17,2]=0
E[17,3]=-dl
E[17,4]=0.76
E[17,5]=0
E[17,6]=17
E[17,7]=(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[17,8]=5
E[17,9]=1

#Element 18 (Variable temperature soil 0-108 cm)
E[18,1]=101
E[18,2]=-dl
E[18,3]=0
E[18,4]=0
E[18,5]=1.08
E[18,6]=18
E[18,7]=(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[18,8]=5
E[18,9]=1

#Element 19 (Variable temperature soil 108-123 cm)
E[19,1]=102
E[19,2]=0.07-dl
E[19,3]=1.08
E[19,4]=0.07
E[19,5]=1.23
E[19,6]=19
E[19,7]=(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[19,8]=5
E[19,9]=1

#Element 20 (Environmental conditions)
E[20,1]=103
E[20,2]=0.07-dl
E[20,3]=1.23
E[20,4]=0.07
E[20,5]=1.85
E[20,6]=20
E[20,7]=(25)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[20,8]=10
E[20,9]=1

#Element 21 (Adiabatic boundary conditions)
E[21,1]=104
E[21,2]=0.07
E[21,3]=1.85
E[21,4]=0.36
E[21,5]=1.85+dl
E[21,6]=21
E[21,7]=(1)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[21,8]=1
E[21,9]=-1

#Element 22 (Internal conditions)
E[22,1]=105
E[22,2]=0.36
E[22,3]=1.27
E[22,4]=0.36+dl
E[22,5]=1.85
E[22,6]=22
E[22,7]=(8)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[22,8]=22
E[22,9]=1

#Element 23 (Internal conditions)
E[23,1]=106
E[23,2]=0.36+dl
E[23,3]=1.27
E[23,4]=1.40
E[23,5]=1.27+dl
E[23,6]=23
E[23,7]=(8)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[23,8]=22
E[23,9]=1

#Element 24 (Adiabatic boundary conditions)
E[24,1]=107
E[24,2]=1.40
E[24,3]=1
E[24,4]=1.40+dl
E[24,5]=1.27
E[24,6]=24
E[24,7]=(1)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[24,8]=1
E[24,9]=-1

#Element 25 (Constant temperature soil)
E[25,1]=108
E[25,2]=0.76+dl
E[25,3]=1-dl
E[25,4]=1.40
E[25,5]=1
E[25,6]=25
E[25,7]=(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[25,8]=5
E[25,9]=1

#Element 26 (Variable temperature soil 0-100 cm)
E[26,1]=109
E[26,2]=0.76
E[26,3]=0
E[26,4]=0.76+dl
E[26,5]=1
E[26,6]=26
E[26,7]=(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[26,8]=5
E[26,9]=1

#Element 27 (Constant temperature soil)
E[27,1]=110
E[27,2]=0
E[27,3]=1.08
E[27,4]=0.07-dl
E[27,5]=1.08+dl
E[27,6]=27
E[27,7]=(1000)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[27,8]=5
E[27,9]=1


text:
Vector{UInt8}(
"POST / HTTP/1.1\r\n" *
"Host: www.somewhere.com\r\n" *
"Content-Length: 9\r\n" *
"\r\n" *
"Body Body"
)


text:
"""

$(TYPEDEF)

Structure that contains some data required to compute the linked cells.

$(TYPEDFIELDS)

"""
struct Box{FloatVector,IntVector}
  sides::FloatVector
  nc::IntVector
  l::FloatVector
  lcell::Int
end

function Box(lcell::Int, sides::AbstractVector, cutoff::Float64)
  # Compute the number of cells in each dimension
  nc = SVector{3,Int}(max.(1,trunc.(Int,sides/(cutoff/lcell))))
  l = SVector{3,Float64}(sides ./ nc)
  Box(sides,nc,l,lcell)
end


text:
module DistributedFlux

using Flux, MPI
using TensorBoardLogger, Logging
using ProgressMeter: next!, Progress
using Flux.Optimise: Params, @progress, batchmemaybe, StopException, update!

include("functor.jl")
include("orthogonal.jl")
include("padding.jl")
include("layers.jl")

call(f, xs...) = f(xs...)
runall(f) = f
runall(fs::AbstractVector) = x -> foreach(call, fs, x)

bcast!(x) = MPI.Initialized() ? MPI.Bcast!(x, 0, MPI.COMM_WORLD) : x

function allreduce!(x)
    MPI.Initialized() || return x
    x′ = zero(x)
    MPI.Allreduce!(x, x′, MPI.SUM, MPI.COMM_WORLD)
    x .= x′ ./ MPI.Comm_size(MPI.COMM_WORLD)
end

unwrap(x) = hasproperty(x, :data) ? x.data : x

function allreduce!(xs, gs)
    for x in xs
        allreduce!(unwrap(gs[x]))
    end
end

function Flux.train!(loss, ps, data, opt, gradient = Flux.gradient; cb_step = x -> (), cb_epoch = x -> (), logger = TBLogger(), verbose = false, sync = false)
    sync && foreach(bcast! ∘ unwrap, ps)
    ps = Params(ps)
    cb_step = runall(cb_step)
    cb_epoch = runall(cb_epoch)
    l̄ = 0f0
    ndata = length(data)
    prog = Progress(ndata)
    @progress for (n, d) in enumerate(data)
        try
            local l
            gs = gradient(ps) do
                l = loss(batchmemaybe(d)...)
            end
            sync && allreduce!(ps, gs)
            update!(opt, ps, gs)
            l̄ = ((n - 1) * l̄ + l) / n
            if verbose
                prog.desc = "$n/$ndata "
                next!(prog, showvalues = [(:loss, l), (:avgloss, l̄)])
                with_logger(logger) do
                    @info "train" loss=l avgloss=l̄
                end
            end
            cb_step((step = n, loss = l, avgloss = l̄))
        catch ex
            if ex isa StopException
                break
            else
                rethrow(ex)
            end
        end
    end
    cb_epoch((avgloss = l̄,))
    return l̄
end

end


text:
abstract type ControllableLoad <: ElectricLoad end

struct InterruptibleLoad <: ControllableLoad
    name::String
    available::Bool
    bus::Bus
    model::String # [Z, I, P]
    maxactivepower::Float64 # [MW]
    maxreactivepower::Float64 # [MVAr]
    sheddingcost::Float64 # $/MWh
    scalingfactor::TimeSeries.TimeArray
end

InterruptibleLoad(; name = "init", status = true, bus = Bus(), model = "0", maxactivepower = 0, maxreactivepower=0, sheddingcost = 999,
                scalingfactor = TimeSeries.TimeArray(Dates.today(),ones(1))) = InterruptibleLoad(name, status, bus, model, maxactivepower, maxreactivepower, sheddingcost, scalingfactor)


text:
cd("c:\\users\\mgjo\\Box Sync\\Documents\\Research projects\\juliaprogram\\constlab")

using Tensors,PGFPlotsX,LaTeXStrings
using ForwardDiff
using DelimetedFiles

include("materials.jl");

model=2

if model==1
#model=1: linear isotropic elasticity
  nstate=1 #cannot set to 0
  npara  =2
  Emod=1000.0; nu=0.3
  mpara=[Emod,nu]
  tol=1e-3
elseif model==2  #von mises, nonlinear mixed hardening

  nstate=10
  npara  =7
  Emod=200000.0; nu=0.3; Sy=100; Hiso=Emod/2; kappainf=30; Hkin=Emod/2; alphainf=50;
  mpara=[Emod,nu,Sy,Hiso,kappainf,Hkin,alphainf]
  tol=1e-3

elseif model==3  #maxwell viscoelasticity

  nstate=1 #cannot set to 0
  npara  =4
  Emod=200000.0; nu=0.3; kappa=1*10^5; mu=1*10^5
  mpara=[Emod,nu,kappa,mu]
  tol=1e-3

end

stateold=zeros(nstate)

#Define load control
#strainctrlv=[1,0,0,0,0,0,0,0,0]
#strain control, uniaxial stress
unkn=[2,3,4,5,6,7,8,9];
notunkn=[1];
#pure strain control
#unkn=[];
#notunkn=[1,2,3,4,5,6,7,8,9];

nunkn=size(unkn,1);
nnotunkn=size(notunkn,1);
stress_unbalance=zeros(nunkn);

tol_unb=1e-3;

#Define time stepping
nsteps=300;
tend=10.; dt=tend/(nsteps-1);
timehist=0.:dt:tend;

#Strain stepping
ϵ_max =0.005; dϵ=ϵ_max/(nsteps-1)
ϵ_hist=0.:dϵ:ϵ_max;

#Stress stepping
σ_max =500.; dσ=(σ_max/(nsteps-1))
σ_hist=0.:dσ:σ_max;
σ_hist=σ_hist.*0.

aa=zeros(nsteps,13)
for i=1:nsteps
   aa[i,:]=[timehist[i],ϵ_hist[i],0.,0.,0.,0.,0.,σ_hist[i],0.,0.,0.,0.,0.]
end
writedlm("loadhist.txt",aa)

ϵ_out = zeros(nsteps) #Float64[]
σ_out = zeros(nsteps) #Float64[]
σ_out2 = zeros(nsteps) #Float64[]

#
ϵ_out,σ_out,σ_out2=time_hist(nsteps,tol_unb,model,tol,σ_hist,ϵ_hist,nstate,stateold,dt,npara,mpara,unkn,notunkn)

println("stress=",σ_out)

#μ = Emod

#=@pgf Axis({thick, no_marks, legend_pos={north_west}, xmajorgrids, ymajorgrids, ylabel=L"\bullet / E", xlabel=L"\epsilon_{11}",width="12cm", height="8cm"},
    Plot(Coordinates(ϵ_hist, σ_out./Emod)),
    LegendEntry(L"\sigma_{11}"),
    Plot(Coordinates(ϵ_out, σ_out2./Emod)),
    LegendEntry(L"\sigma_{22}")
)
=#

@pgf Axis({thick, no_marks, legend_pos={north_west}, xmajorgrids, ymajorgrids, ylabel=L"{\rm [MPa]}", xlabel=L"\epsilon_{11}",width="12cm", height="8cm"},
    Plot(Coordinates(ϵ_hist, σ_out)),
    LegendEntry(L"\sigma_{11}"),
    Plot(Coordinates(ϵ_out, σ_out2)),
    LegendEntry(L"\sigma_{22}")
)


text:
module AlphaGo
using CuArrays
using Flux
using Flux: crossentropy, back!, mse, @treelike
using StatsBase: Weights
using Printf: @sprintf

struct IllegalMove <:Exception end
abstract type Position end
abstract type GameEnv end

struct PositionWithContext
  position::Position
  next_move
  result::Int
end

export MCTSPlayer, NeuralNet, pick_move, play_move!,
      initialize_game!, set_result!,
      should_resign, is_done, result, train!, selfplay,
      train, play, load_model, GoEnv, GomokuEnv, Position, all_legal_moves,
      ChessEnv, result

include("game/go/go.jl")
include("game/gomoku/gomoku.jl")
include("game/env.jl")
include("mcts.jl")
include("mcts_play.jl")
include("features.jl")
include("neural_net.jl")
include("selfplay.jl")
include("train.jl")
include("play.jl")

end #module


text:
module FlatZinc

import MathOptInterface
import ConstraintProgrammingExtensions

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const CleverDicts = MOIU.CleverDicts
const CP = ConstraintProgrammingExtensions

# Formal grammar: https://www.minizinc.org/doc-2.5.5/en/fzn-spec.html#grammar
include("model.jl")
include("export.jl")
include("import.jl")
include("optimizer.jl")

end


text:
__precompile__()

module Postgres

    using DataFrames
    export  PostgresServer,
            PostgresError,
            PostgresServerError,
            PostgresConnection,
            PostgresCursor,
            BufferedPostgresCursor,
            StreamedPostgresCursor,
            PostgresResult,

            query,
            status,
            cursor,
            execute,
            fetch,
            begin_!,
            commit!,
            rollback!,
            escape_value,
            copyto

    #####   interface stub
    # things the developer does not want you to do
    # e.g. call select on a closed connection
    abstract DatabaseError <: Exception
    abstract DatabaseInterface
    abstract DatabaseConnection

    abstract DatabaseCursor
    abstract BufferedDatabaseCursor <: DatabaseCursor
    abstract StreamedDatabaseCursor <: DatabaseCursor

    module Libpq
        include("libpq.jl")
        export PGconn, PGresult, ResultPtr, ConnPtr
    end

    module Types
        include("types.jl")

        export  base_types,
                AbstractPostgresType,
                PostgresType,
                PostgresDomainType,
                PostgresEnumType,
                PostgresValue,
                unsafe_parse,
                naval
    end

    module Results
        using DataArrays
        abstract DatabaseResult
        abstract DatabaseResultInfo
        # things the database does not like e.g. a syntax error in the SQL
        abstract DatabaseServerError <: Exception
        using ..Libpq
        using ..Types
        include("results.jl")

        export PostgresResult,
            PostgresServerError,
            column,
            row
    end

    using .Results
    using .Types
    using .Libpq
    include("postgres.jl")
end

max_stars_repo_name: NCarson/Postgres

max_issues_repo_name: NCarson/Postgres

max_forks_repo_name: NCarson/Postgres


text:
@testset "Test handlers" begin

using HDiscontinuousGalerkin
using Tensors

mesh = rectangle_mesh(TriangleCell, (2,2), Vec{2}((0.0,0.0)), Vec{2}((1.0,1.0)))
dim = 2
Wh = ScalarFunctionSpace(mesh, ContinuousLagrange{dim,RefTetrahedron,1}(); update_face_data = false)
u_h = TrialFunction(Wh)

dh = DofHandler([u_h],mesh)

@test dh.cell_dofs == [1, 2, 3, 2, 4, 3, 2, 5, 4, 5, 6, 4, 3, 4, 7, 4, 8, 7, 4, 6, 8, 6, 9, 8]
@test dh.cell_dofs_offset == [1, 4, 7, 10, 13, 16, 19, 22, 25]
@test dof_range(dh, u_h) == 1:3

# ### Boundary conditions
dbc = Dirichlet(u_h, dh, "boundary", x -> 0)
@test dbc.prescribed_dofs == [1,2,3,5,6,7,8,9]

end


text:
# File for generating data for comparison table

using Pkg
using Random, Distributions
using Statistics
using ProgressMeter
using DataFrames
using JLD2

include("/home/christian/UniMA/EMMA/src/RFR.jl")
include("/home/christian/UniMA/EMMA/src/cross_val.jl")
include("/home/christian/UniMA/EMMA/src/aux_functions.jl")

function big_comp(cv_list, n_list, func="friedman", res=nothing)

    if res===nothing
        # Initiate output DataFrame
        res = DataFrame()
    end

    approaches = ["CART", "weighted", "max depth", "min samples leaf"]

    for n1 in n_list
        x_train, x_test, y_train, y_test = make_data(n1, d, func, σ)

        for (ind, cv) in enumerate(cv_list)
            fit!(cv, x_train, y_train)
            rf = best_model(cv, "mse")
            pred = predict(rf, x_test)
            mse_arr = get_mse(pred, y_test)

            res = vcat(res, DataFrame(bias=mse_arr[1],
                                        variance=mse_arr[2],
                                        mse=mse_arr[3],
                                        Approach=approaches[ind],
                                        n=n1))

        end
    end

    return res
end

### Make comparison
d = 10
σ = 1
m_features = d
cv_parameter = collect(1:5:31)

# Making dictionaries for different approaches
d_base = Dict{Symbol, Vector{Float64}}(
    :max_features => [m_features],
    :n_trees => [30])

d0 = copy(d_base)
d0[:α] = [0.0]

d1 = copy(d_base)
d1[:α] = cv_parameter

d2 = copy(d_base)
d2[:max_depth] = cv_parameter[2:end] # Have to start one above, otherwise depth approach does not work

d3 = copy(d_base)
d3[:min_samples_leaf] = cv_parameter

# Initializing cross validation object for different approaches
cv0 = cross_val(d0, random_state=0)
cv1 = cross_val(d1, random_state=0)
cv2 = cross_val(d2, random_state=0)
cv3 = cross_val(d3, random_state=0)

# List of observations
n_list = [250, 500, 1000, 2000, 4000, 8000, 16000]
# List of CV objects
cv_list = [cv0, cv1, cv2, cv3]
# Repetitions
reps = 100

## Make Friedman simulation
res_friedman = big_comp(cv_list, n_list, "friedman")

for i in 1:reps
    res_friedman = big_comp(cv_list, n_list, "friedman", res_friedman)
    println("\n Done with round", i)
end

jldsave("data/comp_friedman.jld2"; res)

## Make DP3 simulation
res_dp3 = big_comp(cv_list, n_list, "dp3")

for i in 1:reps
    res_dp3 = big_comp(cv_list, n_list, "dp3", res_dp3)
    println("\n Done with round", i)
end

jldsave("data/comp_dp3.jld2"; res_dp3)

## Make DP8 simulation
res_dp8 = big_comp(cv_list, n_list, "dp8")

for i in 1:reps
    res_dp8 = big_comp(cv_list, n_list, "dp8", res_dp8)
    println("\n Done with round", i)
end

jldsave("data/comp_dp8.jld2"; res)

## Make robot simulation
res_robot = big_comp(cv_list, n_list, "robot")

for i in 1:reps
    res_robot = big_comp(cv_list, n_list, "robot", res_robot)
    println("\n Done with round", i)
end

jldsave("data/comp_robot.jld2"; res)

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
#Copyright (c) 2015: Julia Computing Inc. All rights reserved.

using JuliaInXL
using Base.Test

@test xldate(1; use1904windowing=true,roundtoSeconds=true) == Date(1904, 1, 2)
@test xldate(1) == Date(1900, 1,1)
@test xldate(42371.49) == DateTime("2016-01-02T11:45:36")


text:
"""

julia dcgan.jl --outdir ~/dcgan-out
julia dcgan.jl -h # to see all other script options

This example implements a DCGAN (Deep Convolutional Generative Adversarial Network) on MNIST dataset. This implemented model is not identical with the original model. LeNet is used as a base to adapt DCGAN to original MNIST data.1

* Paper url: https://arxiv.org/abs/1511.06434

"""
module DCGAN
using Knet,CUDA,MLDatasets,ArgParse,Printf

function main(args)
    o = parse_options(args)
    o[:seed] > 0 && Knet.seed!(o[:seed])

    # load models, data, optimizers
    wd, wg, md, mg = load_weights(o[:atype], o[:zdim], o[:loadfile])
    xtrn,ytrn = MNIST.traindata()
    xtst,ytst = MNIST.testdata()
    dtrn = minibatch(xtrn, ytrn, o[:batchsize]; shuffle=true, xtype=o[:atype], xsize=(size(xtrn,1),size(xtrn,2),1,o[:batchsize]))
    optd = map(wi->eval(Meta.parse(o[:optim])), wd)
    optg = map(wi->eval(Meta.parse(o[:optim])), wg)
    z = sample_noise(o[:atype],o[:zdim],prod(o[:gridsize]))

    if o[:outdir] != nothing && !isdir(o[:outdir])
        mkpath(o[:outdir])
        mkpath(joinpath(o[:outdir],"models"))
        mkpath(joinpath(o[:outdir],"generations"))
    end

    # training
    println("training started..."); flush(stdout)
    for epoch = 1:o[:epochs]
        dlossval = glossval = 0
        for (x,y) in progress(dtrn)
            noise = sample_noise(o[:atype],o[:zdim],length(y))
            dlossval += train_discriminator!(wd,wg,md,mg,2x .- 1,y,noise,optd,o)
            noise = sample_noise(o[:atype],o[:zdim],length(y))
            glossval += train_generator!(wg,wd,mg,md,noise,y,optg,o)
        end
        dlossval /= length(dtrn); glossval /= length(dtrn)
        println((:epoch,epoch,:dloss,dlossval,:gloss,glossval))
        flush(stdout)

        # save models and generations
        if o[:outdir] != nothing
            filename = @sprintf("%04d.png",epoch)
            filepath = joinpath(o[:outdir],"generations",filename)
            plot_generations(
                wg, mg; z=z, savefile=filepath,
                scale=o[:gridscale], gridsize=o[:gridsize])

            filename = @sprintf("%04d.jld2",epoch)
            filepath = joinpath(o[:outdir],"models",filename)
            save_weights(filepath,wd,wg,md,mg)
        end
    end

    return wd,wg,md,mg
end

function parse_options(args)
    s = ArgParseSettings()
    s.description =
        "Deep Convolutional Generative Adversarial Networks on MNIST."

    @add_arg_table! s begin
        ("--atype"; default="$(Knet.array_type[])"; help="array and float type to use")
        ("--batchsize"; arg_type=Int; default=100; help="batch size")
        ("--zdim"; arg_type=Int; default=100; help="noise dimension")
        ("--epochs"; arg_type=Int; default=20; help="# of training epochs")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
        ("--gridsize"; arg_type=Int; nargs=2; default=[8,8])
        ("--gridscale"; arg_type=Float64; default=2.0)
        ("--optim"; default="Adam(;lr=0.0002, beta1=0.5)")
        ("--loadfile"; default=nothing; help="file to load trained models")
        ("--outdir"; default=nothing; help="output dir for models/generations")
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    o[:atype] = eval(Meta.parse(o[:atype]))
    if o[:outdir] != nothing
        o[:outdir] = abspath(o[:outdir])
    end
    return o
end

function load_weights(atype,zdim,loadfile=nothing)
    if loadfile == nothing
        wd, md = initwd(atype)
        wg, mg = initwg(atype,zdim)
    else
        Knet.@load loadfile wd wg md mg
        wd = convert_weights(wd, atype)
        wg = convert_weights(wg, atype)
        md = convert_moments(md, atype)
        mg = convert_moments(mg, atype)
    end
    return wd, wg, md, mg
end

function save_weights(savefile,wd,wg,md,mg)
    Knet.save(savefile,
              "wd", convert_weights(wd),
              "wg", convert_weights(wg),
              "md", convert_moments(md),
              "mg", convert_moments(mg))
end

function convert_weights(w, atype=Array{Float32})
    w0 = map(wi->convert(atype, wi), w)
    w1 = convert(Array{Any}, w0)
end

function convert_moments(moments,atype=Array{Float32})
    clone = map(mi->bnmoments(), moments)
    for k = 1:length(clone)
        if moments[k].mean != nothing
            clone[k].mean = convert(atype, moments[k].mean)
        end

        if moments[k].var != nothing
            clone[k].var = convert(atype, moments[k].var)
        end
    end
    return convert(Array{Any,1}, clone)
end

function leaky_relu(x, alpha=0.2)
    pos = max.(0,x)
    neg = min.(0,x) .* eltype(x)(alpha)
    return pos .+ neg
end

function sample_noise(atype,zdim,nsamples,mu=0.5,sigma=0.5)
    noise = convert(atype, randn(zdim,nsamples))
    mu = eltype(noise)(mu)
    sigma = eltype(noise)(sigma)
    normalized = (noise .- mu) ./ sigma
end

function initwd(atype, winit=0.01)
    w = Any[]
    m = Any[]

    push!(w, winit*randn(5,5,1,20))
    push!(w, bnparams(20))
    push!(m, bnmoments())

    push!(w, winit*randn(5,5,20,50))
    push!(w, bnparams(50))
    push!(m, bnmoments())

    push!(w, winit*randn(500,800))
    push!(w, bnparams(500))
    push!(m, bnmoments())

    push!(w, winit*randn(2,500))
    push!(w, zeros(2,1))
    return convert_weights(w,atype), m
end

function dnet(w,x0,m; training=true, alpha=0.2)
    x1 = dlayer1(x0, w[1:2], m[1]; training=training)
    x2 = dlayer1(x1, w[3:4], m[2]; training=training)
    x3 = reshape(x2, 800,size(x2,4))
    x4 = dlayer2(x3, w[5:6], m[3]; training=training)
    x5 = w[end-1] * x4 .+ w[end]
end

function dlayer1(x0, w, m; stride=1, padding=0, alpha=0.2, training=true)
    x = conv4(w[1], x0; stride=stride, padding=padding)
    x = batchnorm(x, m, w[2]; training=training)
    x = leaky_relu(x,alpha)
    x = pool(x) #TODO: add mode=2 after it is supported by cuarrays.
end

function dlayer2(x, w, m; training=true, alpha=0.2)
    x = w[1] * x
    x = batchnorm(x, m, w[2]; training=training)
    x = leaky_relu(x, alpha)
end

function dloss(w,m,real_images,real_labels,fake_images,fake_labels)
    yreal = dnet(w,real_images,m)
    real_loss = nll(yreal, real_labels)
    yfake = dnet(w,fake_images,m)
    fake_loss = nll(yfake, fake_labels)
    return real_loss + fake_loss
end

dlossgradient = gradloss(dloss)

function train_discriminator!(wd,wg,md,mg,real_images,ygold,noise,optd,o)
    fake_images = gnet(wg,noise,mg; training=true)
    nsamples = div(length(real_images),784)
    real_labels = ones(Int64, 1, nsamples)
    fake_labels = 2ones(Int64, 1, nsamples)
    gradients, lossval = dlossgradient(
        wd,md,real_images,real_labels,fake_images,fake_labels)
    update!(wd, gradients, optd)
    return lossval
end

function initwg(atype=Array{Float32}, zdim=100, winit=0.01)
    w = Any[]
    m = Any[]

    # 2 dense layers combined with batch normalization layers
    push!(w, winit*randn(500,zdim))
    push!(w, bnparams(500))
    push!(m, bnmoments())

    push!(w, winit*randn(800,500)) # reshape 4x4x16
    push!(w, bnparams(800))
    push!(m, bnmoments())

    # 3 deconv layers combined with batch normalization layers
    push!(w, winit*randn(2,2,50,50))
    push!(w, bnparams(50))
    push!(m, bnmoments())

    push!(w, winit*randn(5,5,20,50))
    push!(w, bnparams(20))
    push!(m, bnmoments())

    push!(w, winit*randn(2,2,20,20))
    push!(w, bnparams(20))
    push!(m, bnmoments())

    # final deconvolution layer
    push!(w, winit*randn(5,5,1,20))
    push!(w, winit*randn(1,1,1,1))
    return convert_weights(w,atype), m
end

function gnet(wg,z,m; training=true)
    x1 = glayer1(z, wg[1:2], m[1]; training=training)
    x2 = glayer1(x1, wg[3:4], m[2]; training=training)
    x3 = reshape(x2, 4,4,50,size(x2,2))
    x4 = glayer2(x3, wg[5:6], m[3]; training=training)
    x5 = glayer3(x4, wg[7:8], m[4]; training=training)
    x6 = glayer2(x5, wg[9:10], m[5]; training=training)
    x7 = tanh.(deconv4(wg[end-1], x6) .+ wg[end])
end

function glayer1(x0, w, m; training=true)
    x = w[1] * x0
    x = batchnorm(x, m, w[2]; training=training)
    x = relu.(x)
end

function glayer2(x0, w, m; training=true)
    x = deconv4(w[1], x0; stride=2)
    x = batchnorm(x, m, w[2]; training=training)
end

function glayer3(x0, w, m; training=true)
    x = deconv4(w[1], x0)
    x = batchnorm(x, m, w[2]; training=training)
    x = relu.(x)
end

function gloss(wg,wd,mg,md,noise,ygold)
    fake_images = gnet(wg,noise,mg)
    ypred = dnet(wd,fake_images,md)
    return nll(ypred, ygold)
end

glossgradient = gradloss(gloss)

function train_generator!(wg,wd,mg,md,noise,labels,optg,o)
    ygold = ones(Int64, 1, length(labels))
    gradients, lossval = glossgradient(wg,wd,mg,md,noise,ygold)
    update!(wg,gradients,optg)
    return lossval
end

function plot_generations(
    wg, mg; z=nothing, gridsize=(8,8), scale=1.0, savefile=nothing)
    if z == nothing
        nimg = prod(gridsize)
        zdim = size(wg[1],2)
        atype = typeof(wg[1]) # wg[1] isa KnetArray ? KnetArray{Float32} : Array{Float32}
        z = sample_noise(atype,zdim,nimg)
    end
    output = Array(0.5 .* (1 .+ gnet(wg,z,mg; training=false)))
    images = map(i->output[:,:,:,i], 1:size(output,4))
    grid = make_image_grid(images; gridsize=gridsize, scale=scale)
    if savefile == nothing
        display(colorview(Gray, grid))
    else
        Knet.save(savefile, grid)
    end
end

splitdir(PROGRAM_FILE)[end] == "dcgan.jl" && main(ARGS)

end # module


text:
# Pretty hacky way of evaluating in this module
struct Engine end

function Base.getproperty(::Engine, s::Symbol)
    prop = py"getattr(eng, $(string(s)))"
    if s == :workspace
        return Workspace(prop)
    else
        return prop
    end
end

Base.propertynames(me::Engine) = Symbol.(py"dir(eng)")

# This is only here to convert outputs of the engine workspace
struct Workspace
    pyobj::PyCall.PyObject
end

Base.getproperty(::Workspace, s::Symbol) = Base.eval(Main, _mat_str(string(s)))

# Not sure why I have to do it this way here but not in the mat" string macro
function Base.setproperty!(::Workspace, s::Symbol, x)
    x_str = _string_for_matlab(x)
    str = string(s) * " = " * x_str
    _mat_str(str)
    return nothing
end

Base.getindex(ws::Workspace, s) = getproperty(ws, Symbol(s))

Base.keys(ws::Workspace) = keys(getfield(ws, :pyobj))

Base.propertynames(ws::Workspace) = propertynames(getfield(ws, :pyobj))

Base.show(io::IO, ws::Workspace) = print(io, getfield(ws, :pyobj))

"""
    connect_matlab(; background=false)
    connect_matlab(engine_name; background=false)

Connect to a named MATLAB session. This must be done before using the integrated MATLAB REPL.
If no name engine name is given, it will connect to the first named session returned from
`find_matlab()`.

## Examples
```julia
julia> using SharedMATLABEngine

julia> connect_matlab(:default_engine);
REPL mode MATLAB initialized. Press > to enter and backspace to exit.

>> a = magic(3)

a =

     8     1     6
     3     5     7
     4     9     2

>> \$a21 = a(2,1)
3.0

julia> a21 .+ mat"a"
3×3 Array{Float64,2}:
 11.0   4.0   9.0
  6.0   8.0  10.0
  7.0  12.0   5.0
```
"""
function connect_matlab(name=nothing; background=false)
    # Note: Don't fix the indentation here, it breaks Julia syntax highlighting
py"""
eng = matlab.engine.connect_matlab($name, background=$background)
"""
    start_repl()
    return Engine()
end

"""
    start_matlab(; option="-nodesktop", background=false)

Start new MATLAB session
"""
function start_matlab(; option="-nodesktop", background=false)
py"""
eng = matlab.engine.start_matlab(option=$option, background=$(background))
"""
    start_repl()
    return Engine()
end

"""
    find_matlab()

Find available MATLAB sessions
"""
find_matlab() = py"matlab.engine.find_matlab()"

function start_repl()
    if isinteractive()
        initrepl(str -> Meta.parse("SharedMATLABEngine.matrepl\"$str\"");
            prompt_text = ">> ",
            start_key = ">",
            prompt_color = :default,
            mode_name = "MATLAB",
        )
    end
    return nothing
end


text:
using TrajectoryOptimization
using Test, SafeTestsets

@time begin
    @time @safetestset "Path Functions" begin include("pathFunctionTests.jl") end
end


text:
function solve_modelARCDiag(H, g, nlp_stop, PData::PDataFact, α::T) where {T}
    # Solve the ARC subproblem once diagonalized into Δ
    # printstyled("On est dans solve_modelARCDiag ⇊ \n", color = :red)
    # Setup the problem
    M = ones(T, size(PData.Δ))

    ϵ1 = sqrt(eps(T)) ./ T(100.0)
    ϵ = T.(ϵ1 .* (1.0 .+ PData.λ))
    λ = max(ϵ, PData.λ .+ ϵ)

    d̃ = -(PData.Δ .+ λ .* M) .\ PData.g̃
    normd̃ = sqrt(d̃ ⋅ d̃)
    seuil = normd̃ ./ α

    # Solve the subproblem (Δ + λ I) d̃ = -g̃ such that λ = ||d̃||/α
    d̃, λ = solve_diag(λ, PData.Δ, PData.g̃, seuil, α, ϵ)

    PData.λ = λ
    # Transform back d̃ into d

    d = AInv(PData, d̃)

    return d, λ
end


text:
@testset "Utilities" begin

    @testset "`solve!` does not return anything" begin
        x = Variable()
        p = satisfy(x >= 0)
        output = solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
        @test output === nothing
    end

    @testset "Show" begin
        x = Variable()
        @test sprint(show, x) == """
        Variable
        size: (1, 1)
        sign: real
        vexity: affine
        $(Convex.show_id(x))"""
        fix!(x, 1.0)
        @test sprint(show, x) == """
        Variable
        size: (1, 1)
        sign: real
        vexity: constant
        $(Convex.show_id(x))
        value: 1.0"""

        @test sprint(show, 2*x) == """
        * (constant; real)
        ├─ 2
        └─ real variable (fixed) ($(Convex.show_id(x)))"""

        free!(x)
        p = maximize( log(x), x >= 1, x <= 3 )

        @test sprint(show, p) == """
        maximize
        └─ log (concave; real)
           └─ real variable ($(Convex.show_id(x)))
        subject to
        ├─ >= constraint (affine)
        │  ├─ real variable ($(Convex.show_id(x)))
        │  └─ 1
        └─ <= constraint (affine)
           ├─ real variable ($(Convex.show_id(x)))
           └─ 3

        status: `solve!` not called yet"""

        x = ComplexVariable(2,3)
        @test sprint(show, x) == """
        Variable
        size: (2, 3)
        sign: complex
        vexity: affine
        $(Convex.show_id(x))"""

        # test `MAXDEPTH`
        # We construct a binary tree of depth >= 3
        # to make sure it gets truncated appropriately.
        x = Variable(2)
        y = Variable(2)
        level3 = hcat(x,y)
        level2 = hcat(level3, level3)
        root = hcat(level2, level2)
        p = minimize(sum(x), root == root)
        @test sprint(show, p) == """
        minimize
        └─ sum (affine; real)
           └─ 2-element real variable ($(Convex.show_id(x)))
        subject to
        └─ == constraint (affine)
           ├─ hcat (affine; real)
           │  ├─ hcat (affine; real)
           │  │  ├─ …
           │  │  └─ …
           │  └─ hcat (affine; real)
           │     ├─ …
           │     └─ …
           └─ hcat (affine; real)
              ├─ hcat (affine; real)
              │  ├─ …
              │  └─ …
              └─ hcat (affine; real)
                 ├─ …
                 └─ …

        status: `solve!` not called yet"""

        # test `MAXWIDTH`
        x = Variable()
        p = satisfy([ x == i for i = 1:100])
        old_maxwidth = Convex.MAXWIDTH[]
        Convex.MAXWIDTH[] = 2
        @test sprint(show, p) == """
            minimize
            └─ 0
            subject to
            ├─ == constraint (affine)
            │  ├─ real variable ($(Convex.show_id(x)))
            │  └─ 1
            ├─ == constraint (affine)
            │  ├─ real variable ($(Convex.show_id(x)))
            │  └─ 2
            ⋮

            status: `solve!` not called yet"""
        Convex.MAXWIDTH[] = old_maxwidth

        # solved problem
        x = Variable()
        p = satisfy(x >= 0)
        output = solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
        @test sprint(show, p) == """
                minimize
                └─ 0
                subject to
                └─ >= constraint (affine)
                   ├─ real variable ($(Convex.show_id(x)))
                   └─ 0

                termination status: OPTIMAL
                primal status: FEASIBLE_POINT
                dual status: FEASIBLE_POINT"""
    end

    @testset "clearmemory" begin
        @test_deprecated Convex.clearmemory()
    end

    @testset "ConicObj with type $T" for T = [UInt32, UInt64]
        c = ConicObj()
        z = zero(T)
        @test !haskey(c, z)
        c[z] = (1, 1)
        @test c[z] == (1, 1)
        x = T[]
        for (k, v) in c
            push!(x, k)
        end
        @test x == collect(keys(c))
        d = copy(c)
        @test d !== c
    end

    @testset "length and size" begin
        x = Variable(2,3)
        @test length(x) == 6
        @test size(x) == (2, 3)
        @test size(x, 1) == 2
        @test size(x, 2) == 3

        x = Variable(3)
        @test length(x) == 3
        @test size(x) == (3, 1)

        x = Variable()
        @test length(x) == 1
        @test size(x) == (1, 1)
    end

    @testset "lastindex and axes" begin
        x = Variable(2, 3)
        @test axes(x) == (Base.OneTo(2), Base.OneTo(3))
        @test axes(x, 1) == Base.OneTo(2)
        @test lastindex(x) == 6
        @test lastindex(x, 2) == 3
        y = x[:,end]
        @test y isa AbstractExpr
        @test size(y) == (2, 1)
    end

    @testset "Parametric constants" begin
        z = Constant([1.0 0.0im; 0.0 1.0])
        @test z isa Constant{Matrix{Complex{Float64}}}

        # Helper functions
        @test Convex.ispos(1)
        @test Convex.ispos(0)
        @test !Convex.ispos(-1)
        @test Convex.ispos([0,1,0])
        @test !Convex.ispos([0,-1,0])
        @test Convex.isneg(-1)
        @test Convex.isneg(0)
        @test !Convex.isneg(1)
        @test Convex.isneg([0,-1,0])
        @test !Convex.isneg([0,1,0])
        @test Convex._size(3) == (1, 1)
        @test Convex._sign(3) == Positive()
        @test Convex._size([-1,1,1]) == (3, 1)
        @test Convex._sign([-1,1,1]) == NoSign()
        @test Convex._sign([-1,-1,-1]) == Negative()
        @test Convex._size([0 0; 0 0]) == (2, 2)
        @test Convex._sign([0 0; 0 0]) == Positive()
        @test Convex._size(0+1im) == (1, 1)
        @test Convex._sign(0+1im) == ComplexSign()

        @test Convex.imag_conic_form(Constant(1.0)) == [0.0]
        @test Convex.imag_conic_form(Constant([1.0, 2.0])) == [0.0, 0.0]
    end

    @testset "Base.vect" begin
    # Issue #223: ensure we can make vectors of variables
    @test size([Variable(2), Variable(3,4)]) == (2,)
    end

    @testset "Iteration" begin
        x = Variable(2,3)
        s = sum([xi for xi in x])
        x.value = [1 2 3; 4 5 6]
        # evaluate(s) == [21] (which might be wrong? expected 21)
        # but [21][1] === 21[1] === 21
        # so this should pass even after "fixing" that
        @test evaluate(s)[1] == 21

        x = Variable(4)
        @test [xi.inds for xi in x] == [1:1, 2:2, 3:3, 4:4]

        x = Variable(0)
        @test [xi for xi in x] == []
        @test iterate(x) == nothing
    end
    # returns [21]; not sure why
    # context("iteration") do
    #     x = Variable(2,3)
    #     s = sum([xi for xi in x])
    #     x.value = [1 2 3; 4 5 6]
    #     @fact evaluate(s) --> 21
    # end
end


text:
@testset "Test Constructors" begin
    @testset "Kalman Filter Constructors" begin
        x₀ = 1
        F = 2
        H = 3
        P₀ = 5
        Q = 7
        R = 11
        kf = KalmanFilter(x₀, F, H, P₀, Q, R)
        @test typeof(kf) == KalmanFilter{Float64}
        # Promotion
        x₀ = [2;0]
        F = Float32.([0.5 1; 0 1])
        B = Int32.([1;0])
        H = [1 0]
        P₀ = I(2)
        Q = 10*I(2)
        R = 0.7*I(1)
        kf = KalmanFilter(x₀, F, H, B, P₀, Q, R)
        @test typeof(kf.B) == Matrix{Float64}
        @test kf.pdy[] == 1

        x₀ = [0.;2.]
        F = [1. 1;0 1]
        H = [1. 0]
        B = Matrix{Float64}(undef, 2,1)
        B[:,:] = [0; 1.]
        P₀ = [10. 0; 0 10.]
        Q = [1. 0; 0 1.]
        R = fill(1.,1,1)
        offset = [1.;1.]
        kf = KalmanFilter(x₀, F, H, B, P₀, Q, R, offset)
        expected = (x₀, F, H, B, P₀, Q, R, offset)
        reality = (kf.x, kf.F, kf.H, kf.B, kf.P, kf.Q, kf.R, kf.offset)

        @test expected==reality
    end

    @testset "Minimax & Bayesian MMAE Constructors" begin
        x₀ = [0.;2.]
        F = [1. 1;0 1]
        H = [1. 0]
        B = Matrix{Float64}(undef, 2,1)
        B[:,:] = [0; 1.]
        P₀ = [10. 0; 0 10.]
        Q = [1. 0; 0 1.]
        R = fill(1.,1,1)
        offset = [1.;1.]
        kf1 = KalmanFilter(x₀, F, H, B, P₀, Q, R, offset)
        kf2 = KalmanFilter(x₀, -F, H, B, P₀, Q, R, offset)
        mini = MinimaxMMAE([kf1, kf2], () ->  Hypatia.Optimizer())
        @test mini.filterbank[1].F == -mini.filterbank[2].F
        bayesian = BayesianMMAE([kf1, kf2])
        @test bayesian.filterbank[1].pdy[] == 1/2
        @test bayesian.filterbank[1].F == -bayesian.filterbank[2].F
    end
end


text:
function GetColElements(T::Type, buff::String, cols::Vararg{Int}; cmt="#")::Array{Array{T, 1}, 1}
    ncols  = length(cols)
    colmin = minimum(cols)
    colmax = maximum(cols)
    cont = Array{Array{T, 1}, 1}()
    for iline in split(buff, "\n")
        words = split(split(iline, cmt)[1])
        if length(words) < colmin
            continue
        end
        elem  = Array{T, 1}()
        for (i, iword) in enumerate(words)
            if colmax < i
                break
            end
            if i in cols
                if T <: Number
                    try
                        push!(elem, parse(T, iword))
                    catch
                        push!(elem, nothing)
                    end
                elseif T == String
                    push!(elem, iword)
                else
                    push!(elem, nothing)
                end
            end
        end
        if 0 < length(elem) <= ncols
            push!(cont, elem)
        end
    end
    return cont
end

max_stars_repo_name: yoshphys/julia-utility

max_issues_repo_name: yoshphys/julia-utility

max_forks_repo_name: yoshphys/julia-utility


text:
module BalanceLaws

using ..VariableTemplates

export BalanceLaw,
    vars_state,
    number_states,
    init_state_prognostic!,
    init_state_auxiliary!,
    compute_gradient_flux!,
    compute_gradient_argument!,
    transform_post_gradient_laplacian!,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_state!,
    update_auxiliary_state!,
    update_auxiliary_state_gradient!,
    nodal_init_state_auxiliary!,
    nodal_update_auxiliary_state!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

include("state_types.jl")
include("interface.jl")

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
#
# Linear fit
#

struct Linear
  a :: Float64
  b :: Float64
  R :: Float64
  x :: Vector{Float64}
  y :: Vector{Float64}
  ypred :: Vector{Float64}
  residues :: Vector{Float64}
end

"""
`fitlinear(x,y)`

Obtains the linear fit: ``y = a*x + b``

Optional lower and upper bounds for a, and constant b can be provided using, for example:

```
fitlinear(x,y, l=lower(a=0.), b=3)
```

# Examples
```jldoctest
julia> x = sort(rand(10)) ; y = sort(rand(10));

julia> fit = fitlinear(x,y)

------------------- Linear Fit -------------

Equation: y = ax + b

With: a = 1.0448783208110997
      b = 0.18817627115683894

Pearson correlation coefficient, R = 0.8818586822210751
Average absolute residue = 0.14274752107157443

Predicted Y: ypred = [0.1987357699444139, 0.32264343301109627...
residues = [0.1613987313816987, 0.22309410865095275...

--------------------------------------------

```
"""
function fitlinear(X :: AbstractArray{<:Real}, Y :: AbstractArray{<:Real};
                   l :: lower = lower(), u :: upper = upper(), b = nothing,
                   options :: Options = Options() )
  # Check data
  X, Y = checkdata(X,Y,options)
  # Set bounds
  vars = [ VarType(:a,Number,1),
           VarType(:b,Nothing,1) ]
  lower, upper = setbounds(vars,l,u)
  if b == nothing
    # Set model
    @. model(x,p) = p[1]*x + p[2]
    # Initial point
    p0 = Vector{Float64}(undef,2)
    initP!(p0,options,lower,upper)
    # Fit
    fit = curve_fit(model, X, Y, p0, lower=lower, upper=upper)
    # Analyze results and return
    R = pearson(X,Y,model,fit)
    x, y, ypred = finexy(X,length(X),model,fit)
    return Linear(fit.param...,R,x,y,ypred,fit.resid)
  else
    lower = [lower[1]]; upper = [upper[1]];
    # Set model
    @. model_const(x,p) = p[1]*x + b
    # Initial point
    p0 = Vector{Float64}(undef,1)
    initP!(p0,options,lower,upper)
    # Fit
    fit = curve_fit(model_const, X, Y, p0, lower=lower, upper=upper)
    # Analyze results and return
    R = pearson(X,Y,model_const,fit)
    x, y, ypred = finexy(X,length(X),model_const,fit)
    return Linear(fit.param...,b,R,x,y,ypred,fit.resid)
  end
end

function Base.show( io :: IO, fit :: Linear )
  println("")
  println(" ------------------- Linear Fit ------------- ")
  println("")
  println(" Equation: y = ax + b ")
  println("")
  println(" With: a = ", fit.a)
  println("       b = ", fit.b)
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println(" Average square residue = ", mean(fit.residues.^2))
  println("")
  println(" Predicted Y: ypred = [",fit.ypred[1],", ",fit.ypred[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" -------------------------------------------- ")
end

export fitlinear


text:
"""
    ReachabilitySolver{S} <: Solver

Solves reachability and constrained reachability problems in MDPs and POMDPs.
It returns the policy that maximizes the probability of reaching a given set of states.
It takes as input the set of states to reach and the set of states to avoid, as well as the underlying solver.
Any solver from POMDPs.jl are supported.

# Fields

The field are specified as keyword arguments to the solver.

- `reach::Set{S}` the set of states to reach
- `avoid::Set{S}` the set of states to avoid
- `solver::Solver` the underlying solver to use (default is `ValueIterationSolver`)

"""
@with_kw struct ReachabilitySolver{S} <: Solver
    reach::Set{S} = Set{S}()
    avoid::Set{S} = Set{S}()
    solver::Solver = ValueIterationSolver()# can use any solver that returns a value function :o
end

function POMDPs.solve(solver::ReachabilitySolver, mdp::MDP)
    rmdp = ReachabilityMDP(mdp, solver.reach, solver.avoid)
    policy = solve(solver.solver, rmdp)
    return ReachabilityPolicy(policy, mdp, solver.reach, solver.avoid)
end

function POMDPs.solve(solver::ReachabilitySolver, pomdp::POMDP)
    rmdp = ReachabilityPOMDP(pomdp, solver.reach, solver.avoid)
    policy = solve(solver.solver, rmdp)
    return ReachabilityPolicy(policy, pomdp, solver.reach, solver.avoid)
end

struct ReachabilityPolicy{P <: Policy, M <: Union{MDP, POMDP}, S} <: Policy
    policy::P
    problem::M
    reach::Set{S}
    avoid::Set{S}
end

POMDPs.action(policy::ReachabilityPolicy, s) = action(policy.policy, s)
function POMDPs.value(policy::ReachabilityPolicy, s)
    if s ∈ policy.reach
        return 1.0
    else
        return value(policy.policy, s)
    end
end
function POMDPPolicies.actionvalues(policy::ReachabilityPolicy, s)
    if s ∈ policy.reach
        return ones(length(actions(policy.problem)))
    else
        return actionvalues(policy.policy, s)
    end
end

struct ReachabilityPOMDP{S, A, O, P <: POMDP{S,A,O}} <: POMDP{S,A,O}
    problem::P
    reach::Set{S}
    avoid::Set{S}
end

struct ReachabilityMDP{S, A, P <: MDP{S,A}} <: MDP{S,A}
    problem::P
    reach::Set{S}
    avoid::Set{S}
end

POMDPs.states(r::Union{ReachabilityMDP, ReachabilityPOMDP}) = states(r.problem)
POMDPs.stateindex(r::Union{ReachabilityMDP, ReachabilityPOMDP}, s) = stateindex(r.problem, s)
POMDPs.statetype(r::Union{ReachabilityMDP, ReachabilityPOMDP}) = statetype(r.problem)
POMDPs.actions(r::Union{ReachabilityMDP, ReachabilityPOMDP}) = actions(r.problem)
POMDPs.actionindex(r::Union{ReachabilityMDP, ReachabilityPOMDP}, a) = actionindex(r.problem, a)
POMDPs.actiontype(r::Union{ReachabilityMDP, ReachabilityPOMDP}) = actiontype(r.problem)
POMDPs.observation(r::ReachabilityPOMDP, s, a) = observation(r.problem, s, a)
POMDPs.observations(r::ReachabilityPOMDP) = observations(r.problem)
POMDPs.obsindex(r::ReachabilityPOMDP, o) = obsindex(r.problem, o)
POMDPs.obstype(r::ReachabilityPOMDP) = obstype(r.problem)

POMDPs.transition(r::Union{ReachabilityMDP, ReachabilityPOMDP}, s, a) = transition(r.problem, s, a)

POMDPs.isterminal(r::Union{ReachabilityMDP, ReachabilityPOMDP}, s) = (s ∈ r.reach) || (s ∈ r.avoid)
POMDPs.initialstate_distribution(r::Union{ReachabilityMDP, ReachabilityPOMDP}) = initialstate_distribution(r.problem)
POMDPs.initialstate(r::Union{ReachabilityMDP, ReachabilityPOMDP}, rng::AbstractRNG) = initialstate(r.problem, rng)

POMDPs.reward(r::Union{ReachabilityMDP, ReachabilityPOMDP}, s) = float(s ∈ r.reach)
POMDPs.reward(r::Union{ReachabilityMDP, ReachabilityPOMDP}, s, a) = reward(r, s)
POMDPs.reward(r::Union{ReachabilityMDP, ReachabilityPOMDP}, s, a, sp) = reward(r, sp)
POMDPs.discount(r::Union{ReachabilityMDP, ReachabilityPOMDP}) = 1.0


text:
import NodeJS
npm = NodeJS.npm_cmd()
node_root = normpath(joinpath(@__DIR__, "../node"))
println(run(`$npm --prefix=$node_root --scripts-prepend-node-path=true install $node_root`))

max_stars_repo_name: JuliaPluto/pluto-export

max_issues_repo_name: JuliaPluto/pluto-export

max_forks_repo_name: JuliaPluto/pluto-export


text:
using JLD
using Plots
using LaTeXStrings
using Formatting
using TexTables
#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.75)
#Extracting stored varibles

lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

ld = load("./out/jld/FFTVaryRho_lambda_10.jld")
for key in keys(ld)
  symb = Symbol(key)
  data = ld[key]
  @eval $(symb) = $(data)
end
accplot = plot(xlabel=L"\rho", ylabel=L"\log_{10}\left(\mathrm{RMSE}\right)", legend=:bottomright, ylims=(-4.25,-0.75), xlims=(0.5, 10.0), xticks=1:1:10)

covplot = plot(xlabel=L"\rho", ylabel=L"\mathrm{coverage}", legend=:bottomright, ylims=(0.72, 0.92), xlims=(0.5, 10.0), xticks=1:1:10)

hline!(covplot, [0.9], color=:black, label=L"\mathrm{theoretical}\  \mathrm{coverage}")

plot!(accplot,
      ρList,
      log10.(RMSE_int),
      label=L"\lambda = 1.0, \mathrm{scattered}",
      markershape=:circle,
      markersize=5,
      color=silver,
      linestyle=:solid,
      linewidth=5)

plot!(accplot,
      ρList,
      log10.(RMSE_exp),
      label=L"\lambda = 1.0, \mathrm{region}",
      markershape=:dtriangle,
      markersize=5,
      color=lightblue,
      linestyle=:solid,
      linewidth=5)

plot!(covplot,
      ρList,
      coverage_int,
      label=L"\lambda = 1.0, \mathrm{scattered}",
      markershape=:circle,
      markersize=5,
      color=silver,
      linestyle=:solid,
      linewidth=5)

plot!(covplot,
      ρList,
      coverage_exp,
      label=L"\lambda = 1.0, \mathrm{region}",
      markershape=:dtriangle,
      markersize=5,
      color=lightblue,
      linestyle=:solid,
      linewidth=5)

ld = load("./out/jld/FFTVaryRho_lambda_13.jld")
for key in keys(ld)
  symb = Symbol(key)
  data = ld[key]
  @eval $(symb) = $(data)
end

plot!(accplot,
      ρList,
      log10.(RMSE_int),
      label=L"\lambda = 1.3, \mathrm{scattered}",
      markershape=:circle,
      markersize=5,
      color=rust,
      linestyle=:dash,
      linewidth=5)

plot!(accplot,
      ρList,
      log10.(RMSE_exp),
      label=L"\lambda = 1.3, \mathrm{region}",
      markershape=:utriangle,
      markersize=5,
      color=orange,
      linestyle=:dot,
      linewidth=5)

plot!(covplot,
      ρList,
      coverage_int,
      label=L"\lambda = 1.3, \mathrm{scattered}",
      markershape=:circle,
      markersize=5,
      color=rust,
      linestyle=:dash,
      linewidth=5)

plot!(covplot,
      ρList,
      coverage_exp,
      label=L"\lambda = 1.3, \mathrm{region}",
      markershape=:utriangle,
      markersize=5,
      color=orange,
      linestyle=:dot,
      linewidth=5)

savefig(accplot, "./out/plots/FFT_accplot.pdf")
savefig(covplot, "./out/plots/FFT_covplot.pdf")

# # ----------------------------------------
# # Creating plot of test set
# # ----------------------------------------
#
# xTestPlot = scatter(xTest[1,:], xTest[2,:],
#                     markersize=0.5,
#                     markerstrokewidth=0.1,
#                     legend=false)
#
# savefig(xTestPlot, "./out/plots/FFTxTestScatter.png")

# # ----------------------------------------
# # Creating table
# # ----------------------------------------
#
# outputTables = Array{IndexedTable{1,1},1}(undef, length(ρList))
# for (ρInd, ρ) in enumerate(ρList)
#   global outputTables
#   global coverage
#   global RMSE
#   RMSETab = TableCol("\$\\rho = $ρ\$",
#                      [Symbol("\$\\mathrm{RMSE}\$")],
#                      [FormattedNumber(RMSE[ρInd])])
#
#   coverageTab = TableCol("\$\\rho = $ρ\$",
#                          [Symbol("\$\\mathrm{Coverage}\$")],
#                          [FormattedNumber(coverage[ρInd])])
#   outputTables[ρInd] = vcat(RMSETab, coverageTab)
# end
#
# tab = hcat(outputTables[:]...)
# f = open("./out/tables/FFTTable.tex", "w")
# write(f, to_tex(tab))
# close(f)


text:
g = Multigraph([0 1 0; 1 0 1; 0 1 0])
ps = [Rational(0) for i = 1:3]
v_t = [SpiderType.X, SpiderType.Z, SpiderType.X]
zxd = ZXDiagram(g, v_t, ps)
zxd2 = ZXDiagram(g, Dict(zip(1:3,v_t)), Dict(zip(1:3,ps)))
@test zxd.mg == zxd2.mg && zxd.st == zxd2.st && zxd.ps == zxd2.ps

zxd2 = copy(zxd)
@test zxd.st == zxd2.st && zxd.ps == zxd2.ps
@test ZXCalculus.spider_type(zxd, 1) == SpiderType.X
@test nv(zxd) == 3 && ne(zxd) == 2

@test rem_edge!(zxd, 2, 3)
@test outneighbors(zxd, 2) == inneighbors(zxd, 2)

ZXCalculus.add_spider!(zxd, SpiderType.H, 0//1, [2, 3])
ZXCalculus.insert_spider!(zxd, 2, 4, SpiderType.H)
@test nv(zxd) == 5 && ne(zxd) == 4

zxd3 = ZXDiagram(3)
ZXCalculus.insert_spider!(zxd3, 1, 2, SpiderType.H)
pushfirst_gate!(zxd3, Val{:SWAP}(), [1, 2])
push_gate!(zxd3, Val{:SWAP}(), [2, 3])
@test ZXCalculus.qubit_loc(zxd3, 1) == ZXCalculus.qubit_loc(zxd3, 2) == ZXCalculus.qubit_loc(zxd3, 7)

@testset "float to rational" begin
    @test ZXCalculus.continued_fraction(2.41, 10) === 241//100
    @test ZXCalculus.continued_fraction(1.3, 10) === 13//10
    @test ZXCalculus.continued_fraction(0, 10) === 0//1
    @test ZXCalculus.continued_fraction(-0.5, 10) === -1//2
    zxd = ZXDiagram(4)
    push_gate!(zxd, Val(:X), 3, 0.5)
    @test zxd.ps[9] == 1//2
    push_gate!(zxd, Val(:X), 3, -0.5)
    @test zxd.ps[10] == 3//2
    push_gate!(zxd, Val(:Z), 3, 0)
    @test zxd.ps[11] == 0//1
    @test_warn "" push_gate!(zxd, Val(:Z), 3, sqrt(2))
    @test_throws MethodError push_gate!(zxd, Val(:Z), 3, sqrt(2); autoconvert=false)
end

zxd4 = ZXDiagram(2)
ZXCalculus.add_global_phase!(zxd4, ZXCalculus.Phase(1//2))
ZXCalculus.add_power!(zxd4, 2)
@test ZXCalculus.scalar(zxd4) == Scalar(2, 1//2)


text:
"""Golden Ratio"""
G = (1 + sqrt(5))/2

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
### Gali's Book Chapter 3 model
### We use this as a consistency check of gensys implementation
### By Daniel Coutinho and Gilberto Boaretto

#### Parameters from Gali's book example, chapter 3, page 52, first edition

bet = 0.99
sig = 1
phi = 1
alfa = 1/3
epsilon = 6
theta = 2/3
phi_pi = 1.5
phi_y = 0.5/4
rho_v = 0.5

true_pars = (alfa = alfa, bet = bet, epsilon = epsilon, theta = theta, sig = sig, s2 = 1, phi = phi, phi_pi = phi_pi, phi_y = phi_y, rho_v = rho_v)

THETA = (1-alfa)/(1-alfa+alfa*epsilon)
lamb = (1-theta)*(1-bet*theta)/theta*THETA
kappa = lamb*(sig+(phi+alfa)/(1-alfa))

## Gali's give analytical expressions for the effects of a shock in the monetary policy over inflation and output gap,chapter 3, page 51, first edition

LAMBDA_v = ((1-bet*rho_v)*(sig*(1-rho_v)+phi_y)+kappa*(phi_pi-rho_v))^(-1)

y_irf(v) = -(1-bet*rho_v)*LAMBDA_v*v
pi_irf(v) = -kappa*LAMBDA_v*v
r_irf(v) = sig*(1-rho_v)*(1-bet*rho_v)*LAMBDA_v*v
i_irf(v) = (sig*(1-rho_v)*(1-bet*rho_v) - rho_v*kappa)*LAMBDA_v*v
v_irf(v,uu) = rho_v*v + uu

## Test drive: lets see if this gives the same IRFs that we get in Gali's book
## This is just for debugging, so comment it out when using it to check gensys consistency

#using Plots

irfs_true = zeros(15,5)

irfs_true[1,1] = 0

uu = 0.25

for t in 2:15
    irfs_true[t,4] = v_irf(irfs_true[t-1,4],uu)
    irfs_true[t,1] = pi_irf(irfs_true[t,4])
    irfs_true[t,2] = y_irf(irfs_true[t,4])
    irfs_true[t,3] = i_irf(irfs_true[t,4])
    irfs_true[t,5] = r_irf(irfs_true[t,4])
    global uu = 0
end

#plot(irfs[2:13,1], label = "v")
#plot(irfs[2:13,2]*4, label = "pi")
#plot(irfs[2:13,3], label = "y")
#plot(4*irfs[2:13,4], label = "i")
#plot(4*irfs[2:13,5], label = "r")

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
module PETSc

using MPI, LinearAlgebra, SparseArrays

MPI.Initialized() || MPI.Init()

using Libdl

function _petsc_link(fname)
"""
[`$fname`](https://petsc.org/release/docs/manualpages/$fname.html)
"""
end

function _doc_external(fname)
"""
- PETSc Manual: $(_petsc_link(fname))
"""
end

include("const.jl")
include("startup.jl")
include("lib.jl")
include("init.jl")
include("ref.jl")
include("viewer.jl")
include("options.jl")
include("vec.jl")
include("mat.jl")
include("matshell.jl")
include("dm.jl")
include("dmda.jl")
include("ksp.jl")
include("pc.jl")
include("snes.jl")
include("sys.jl")

end