text:
using Pkg,LinearAlgebra, SparseArrays, Random,Distributions;

y=[1,2,3]
function get_Y(y)
    m=length(y)#number of rows
    n=m*(m-1) #number of columns
    res = zeros(m,n) #set the result
    start = 1
    for i=1:length(y)
        row = copy(y)
        row = deleteat!(row, i)
        res[i,start:start+m-2]=row
        start = start+m-1
    end
    return res
end
get_Y(y)

#fake data
y1=[1,2,3]
y2=[4,5,6]
Y1=get_Y(y1)
Y2=get_Y(y2)
Y=[Y1,Y2]
R=[1 0 0;0 2 0; 0 0 3]
w1=[6,7,8]
w2=[6,3,4]
W=[w1,w2]
t_square = 3
lambda0 = 4
get_Y(y1)

function get_para(Y,R,t_square,lambda0)
    m = size(Y[1])[1]#number of rows
    n=m*(m-1) #number of columns
    first = zeros(n,n)
    second=zeros(n,1)
    for i=1:length(Y)
        first += Y[i]'*inv(R)*Y[i]
        second += Y[i]'*inv(R)*W[i]
    end
    first += Diagonal(repeat([t_square],n))
    second += repeat([lambda0*t_square],n)
    mu = vec(inv(first)*second)
    return mu,Symmetric(inv(first))
end

#MCMC
mu,var = get_para(Y,R,t_square,lambda0)
rand(MvNormal(mu,Symmetric(var)))

y = [1,2,3]
function get_my_Y(y)
    m=length(y)#number of rows
    n::Int64=(m*(m-1))/2 #number of columns
    res = zeros(m,n) #set the result
    start = 1
    for i = 1:length(y)-1
        row = y[2:i+1]
        res[i+1,start:start + i - 1] = row
        start = start + i
    end
    return res
end
get_my_Y(y)

#Making fake data
y1=[1,2,3]
y2=[4,5,6]
Y1=get_my_Y(y1)
Y2=get_my_Y(y2)
Y=[Y1,Y2]
R=[1 0 0;0 2 0; 0 0 3]
w1=[6,7,8]
w2=[6,3,4]
W=[w1,w2]
t_square = 3
lambda0 = 4

function get_para(Y,R,t_square,lambda0,W)
    m = size(Y[1])[1]#number of rows
    n::Int64=(m*(m-1))/2 #number of columns
    first = zeros(n,n)
    second=zeros(n,1)
    for i=1:length(Y)
        first += Y[i]'*inv(R)*Y[i]
        second += Y[i]'*inv(R)*W[i]
    end
    first += Diagonal(repeat([t_square],n))
    second += repeat([lambda0*t_square],n)
    mu = vec(inv(first)*second)
    return mu,Symmetric(inv(first))
end

@time mu,var = get_para(Y,R,t_square,lambda0,W)
#rand(MvNormal(mu,Symmetric(var)))

mu

var

res = randn(3)
res = L.U*res + mu
function sampling(mu,var,n)
    res = randn(n)
    L = cholesky(var)
    res = L.L*res + mu
    return res
end
sampling(mu,var,3)

#Making fake data
y1=[1,2,3]
y2=[4,5,6]
Y1=get_my_Y(y1)
Y2=get_my_Y(y2)
Y=[Y1,Y2]
R=[1 0 0;0 2 0; 0 0 3]
w1=[6 7 8]
w2=[6 3 4]
W=[w1 w2]
t_square = 3
lambda0 = 4
number=2

#Y is the data for all individuals

#unlist function will tranfrom y from a list to a big matrix
function unlist(Y)
    res = Y[1]
    for i = 2:length(Y)
        res = vcat(res,Y[i])
    end
    return res
end

unlist(Y)

function get_para(Y,R,W,t_square,lambda0,number)
    m = size(Y[1])[1]#number of rows
    n::Int64=(m*(m-1))/2 #number of columns
    Y = unlist(Y) #Set big Y matrix

    #Define sparse matrix
    I = collect(1:number*m); J=collect(1:number*m);V = repeat(diag(R),number)
    R = sparse(I,J,V)

    #formula calculation
    first = Y'*inv(Matrix(R))*Y + Diagonal(ones(m,m))
    second = Y'*inv(Matrix(R))*transpose(vcat(W)) + ones(m,1)*lambda0*t_square
    mu = vec(inv(first)*second)
    return mu,Symmetric(inv(first))
end

a, b = get_para(Y,R,W,t_square,lambda0,number)

#Y is a data for single persion
y1 = [1,2,3]
function get_sparse_Y(Y)
    m = size(Y)[1] #number of rows
    n::Int64=(m*(m-1))/2  #number of columns
    col_index = collect(1:n)
    row_index = [];value = []
    for i = 2:m
        current = fill(i,i-1) #get a element i with i-1 time
        row_index = vcat(row_index,current)
        current_value = Y[2:i]
        value = vcat(value,current_value)
    end

    #Return the sparse matrix
    row_index = convert(Array{Int64,1}, row_index) #change Any type
    value = convert(Array{Float64,1}, value) #Change Any type

    res = sparse(row_index,col_index,value)
end

get_sparse_Y(y1)

#Fake data
y1=[1,2,3]
y2=[4,5,6]
ys1=get_sparse_Y(y1)
ys2=get_sparse_Y(y2)
YS = [ys1,ys2]
R=[1 0 0;0 2 0; 0 0 3]
w1=[6,7,8]
w2=[6,3,4]
W=[w1,w2]
t_square = 3
lambda0 = 4
m = 3
number = 10000

YS = repeat([get_sparse_Y(rand(3))],10000)
YS = unlist(YS)
#W=unlist(W)
W = rand(30000)

function get_sparse_para(YS,R,W,t_square,lambda0,number,m)

    #Define residual matrix
    value = diag(R).^(-1)
    RS = kron(sparse(I,number,number),sparse(1:m,1:m,value))

    #formula calculation
    YS_RS_product = YS'*RS
    first = YS_RS_product*YS + t_square*sparse(I,m,m)
    second = YS_RS_product*W + ones(m,1)*lambda0*t_square

    #Return the result
    first = inv(Matrix(first))
    mu = vec(first*second)
    var = Symmetric(first)

    #Take sample
    res = rand(MvNormal(mu,var))
    return res
end

#a,b = get_sparse_para(YS,R,W,t_square,lambda0,number,3)#
@time lambda= get_sparse_para(YS,R,W,t_square,lambda0,number,m)

y1 = rand(1000);y2 = rand(1000);y3=rand(1000)
BigY = [y1;y2;y3]
Lambda = [1 0 0; lambda[1] 1 0;lambda[2] lambda[3] 1]

function Modify_Y(Lambda,Y,n_individuals)
    Identy = Matrix{Int32}(I,n_individuals,n_individuals)
    res = kron(Identy,Lambda)*BigY
    return res
end

Modify_Y(Lambda,BigY,1000)


text:
include("../src/JOH.jl")
include("../MDMKP/MDMKP.jl")
using CPLEX, JuMP, Plots

problems = MDMKP.load_folder()

function graph_problem(problem)
    m = MDMKP.create_MIPS_model(problem, time_limit=1)
	set_optimizer(m, with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=false))
	JOH.Matheur.set_tolerance!(m, .01)
	JOH.Matheur.set_threads!(m, 5)
    gaps = []
    while true
        optimize!(m)
		gap = MOI.get(m, MOI.RelativeGap())
		push!(gaps, gap)
		println(termination_status(m), gap)
		if termination_status(m) == MOI.OPTIMAL
			return gaps
		end
    end
end

looser_gaps = graph_problem.(problems[1:20])

plot(looser_gaps)
plot(gaps)

ext: jl

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
"""
    @test_formany [ntests = 100] [maxtests = 1000] [logto = ""] [argument_data, proposition]

Takes proposition depending on some variables, generates sample of their values
and tests the proposition. Returns `Pass` if proposition holds for all values
in the sample, `Fail` if proposition does not hold for some value, and `Error`
if a part of code could not be evaluated.

# Arguments

Macro accepts a tuple of expressions that describe either statement of a test
or specify optional keyword arguments.

Statement of a test must be a sequence of comma-separated expressions, where
all but last expressions are read as `argument_data` and last statement is
`proposition` depending on these arguments.

* `argument_data` consists of arbitrarily many expressions containing
information about variables, their types and conditions on them. Expressions
must have one of the following forms:

  * `a<x::T<b` or `a>x::T>b`, where middle variable `x` must be a symbol;
  `T`  must be a type; bounds `a` & `b` can be arbitrary expressions;
  * `a<<x::T<b`, `a<x::T<<b`, `a>>x::T>b` or `a>x::T>>b`,
  similar to the one above;
  * shorthand `x` instead of `x::T` may be used, e.g. `a<x<b`,
  then type of variable defaults to `Float64`.

* `proposition`: is the expression to be tested. Two types of propositions are allowed:

  * `proposition` of form `condition-->statement`, where `condition` and
  `statement` are boolean-valued expressions
  * `proposition` is boolean-valued expression itself

Optional keyword arguments include the following:

* `ntests = 100`: number of tests that must be performed.

* `maxtests = 1000`: maximal number of attempts to generate arguments
satisfying `condition`

* `logto = "path_to_file"`: if provided then parsed as a path to file,
where log of tests will be recorded.

# Description
Macro @test_formany is a test that attempts to generate sample of argument values
of size `ntests` satisfying the `condition` if it is given. Then it checks whether
the proposition holds for all elements in the sample. If number of attempts
exceeds `maxtests` parameter before at least `ntests` values for arguments
are generated, macro gives output `Base.Test.Error`.

If `proposition` is of form `condition-->statement` (see arguments list),
then test checks whether implication holds. Namely, the test passes if whenever
`condition` evaluates to `true` for generated argument values, then `statement`
also evaluates to `true`.

If `proposition` is just a boolean-valued expression, code simply checks whether
`prop` holds for generated values.

Some caution is required to avoid vacuous conditions. If test is unable to
generate values satisfying the `condition` after `maxtests` attempts, output
of the test is Base.Test.Error as mentioned before.

Default for optional keyword arguments is the following: `ntests = 100, maxtests = 1000`.
If only `ntests` parameter is specified, `maxtests = 10*ntests`. If
`logto = "path_to_file"` is given, log of test will be written to corresponding
file,for example:

* `@test_formany ntests = 10 maxtests = 1000  logto =  "./test_log.csv" 1<x::Float32<10, proposition`
* `@test_formany ntests = 50 logto = "./test_log.csv" 1<x<10, proposition` ## maxtests = 500
* `@test_formany logto = "./test_log.csv" 1<x<<Inf, proposition` ## ntests = 100; maxtests = 1000
* `@test_formany ntests = 50 maxtests = 1000 1<x<10, proposition` ## not logged
* `@test_formany 1<x::Int32<Inf, proposition logto = "./test_log.csv" maxtests = 1000 ntests = 10` ## order is not important

# Examples
```julia
julia> @test_formany 100>x::Float64>10, x+10<y::Float64<<1000, y-5<z::Float64<Inf, z>x+5
Test Passed
  Expression: (:((100 > x::Float64 > 10,x + 10 < y::Float64 << 1000,y - 5 < z::Float64 < Inf,z > x + 5)),)

julia> @test_formany ntests = 1000 100>x::Float64>10, x+10<y::Float64<<1000, y-5<z::Float64<Inf, z>x+6
Test Failed
  Expression: (:(ntests = 1000),:((100 > x::Float64 > 10,x + 10 < y::Float64 << 1000,y - 5 < z::Float64 < Inf,z > x + 6)))

## example of bad design for generators:
julia> @test_formany -1000<x::Float64<1000, x>999-->x+1>1000
Error During Test
  Test threw an exception of type ErrorException
  Expression: (:((-1000 < x::Float64 < 1000,$(Expr(:-->, :(x > 999), :(x + 1 > 1000))))),)
  Found only 1/100 values satisfying given condition.

## log is increasing
julia> @test_formany Inf>x::Float64>0,Inf>y::Float64>0, x<y-->log(x)<log(y) ntests = 1000
Test Passed
  Expression: (:((Inf > x::Float64 > 0,Inf > y::Float64 > 0,$(Expr(:-->, :(x < y), :(log(x) < log(y)))))),:(ntests = 1000))

## f(x) = x^3 is not convex on [-100,100]
julia> @test_formany 0<a<1,-100<x<100,-100<y<100,(a*x+(1-a)*y)^3<a*x^3+(1-a)*y^3
Test Failed
  Expression: (:((0 < a < 1,-100 < x < 100,-100 < y < 100,(a * x + (1 - a) * y) ^ 3 < a * x ^ 3 + (1 - a) * y ^ 3)),)

## f(x) = x^3 is convex on [-100,100]
julia> @test_formany 0<a<1,0<x<100,0<y<100, (a*x+(1-a)*y)^3<=a*x^3+(1-a)*y^3
Test Passed
  Expression: (:((0 < a < 1,0 < x < 100,0 < y < 100,(a * x + (1 - a) * y) ^ 3 <= a * x ^ 3 + (1 - a) * y ^ 3)),)

## Cobb-Douglas utility function is concave
julia> function respects_concavity(f,t,x1,x2)
                  f(t*x1+(1-t)*x2) > t * f(x1) + (1-t) * f(x2)
              end

julia> function respects_weak_concavity(f,t,x1,x2)
                  f(t*x1+(1-t)*x2) >= t * f(x1) + (1-t) * f(x2)
              end

julia> function cobb_douglas(a,x)
                  prod(x.^a)
              end

# This passes the @test_formany, since probability of discovering a counter-example with a = .5 is zero.
julia> @test_formany 0<t<1, 0<a<1, 1<x1<10, 1<x2<10, 1<y1<10, 1<y2<10,respects_concavity(x -> cobb_douglas([a,1-a],x),t,[x1,y1],[x2,y2])
Test Passed
  Expression: (:((0 < t < 1,0 < a < 1,1 < x1 < 10,1 < x2 < 10,1 < y1 < 10,1 < y2 < 10,respects_concavity((x->begin  # REPL[5], line 1:
                    cobb_douglas([a,1 - a],x)
                end),t,[x1,y1],[x2,y2]))),)

julia> @test_formany 0<t<1, 0<a<1, 1<x1<10, 1<x2<10, 1<y1<10, 1<y2<10,respects_weak_concavity(x -> cobb_douglas([a,1-a],x),t,[x1,y1],[x2,y2])
Test Passed
  Expression: (:((0 < t < 1,0 < a < 1,1 < x1 < 10,1 < x2 < 10,1 < y1 < 10,1 < y2 < 10,respects_weak_concavity((x->begin  # REPL[6], line 1:
                    cobb_douglas([a,1 - a],x)
                end),t,[x1,y1],[x2,y2]))),)

## handling errors
julia> @test_formany -1<x<1, abs(log(x))>0
Error During Test
  Test threw an exception of type DomainError
  Expression: (:((-1 < x < 1,abs(log(x)) > 0)),)
  DomainError
```
"""
macro test_formany(exprs...)
   inner_ex = try
	maxtests = 0; ntests = 0
	prop = :(); cond = :()
	logging  = ""
	var_data = Array(Any,1,6)
	for ex in exprs
		if isa(ex, Expr) && ex.head == :tuple
			var_data = parse_argument_data(ex)
			if ex.args[end].head == :-->
				prop = ex.args[end].args[2]
				cond = ex.args[end].args[1]
			else
				prop = ex.args[end]
				cond = true
			end
		elseif isa(ex, Expr) && ex.head == :(=)
			if ex.args[1] == :ntests
				ntests = ex.args[2]
			elseif ex.args[1] == :maxtests
				maxtests = ex.args[2]
			elseif ex.args[1]== :logto
				logging = ex.args[2]
			else
				error("Invalid macro input $ex.")
			end
		else
			error("Invalid macro input $ex.")
		end
	end
	if ntests == 0
		ntests = 100
	end
	if maxtests == 0
		maxtests = 10*ntests
	end
	num_of_vars = size(var_data,1)
	generate_values = :()
	values = :([])
	for i in 1:num_of_vars
		next_expr = :($(esc(var_data[i,1]))=
			   custom_generator($(esc(var_data[i,2])),$(esc(var_data[i,3])),
			      $(esc(var_data[i,4])),$(esc(var_data[i,5])),$(esc(var_data[i,6])))(div(n,2)+3))
		generate_values = Expr(:block,generate_values.args...,next_expr)
		values = Expr(:vect, values.args..., :($(esc(var_data[i,1]))))
	end
#####
	if logging == ""    ## no logging
		inner_ex = quote try
			for n in 1:$maxtests
				$generate_values
				if $(esc(cond))
					num_good_args += 1
					res = res && $(esc(prop))
					if !res
						fail_data = $(esc(prop))
						break
					end
					if num_good_args >= $ntests
						break
					end
				end
			end
		        # no counter-example but insufficient amount of tests
			if res && num_good_args<$ntests
				nt = $ntests
 				error("Found only $num_good_args/$nt values satisfying given condition.")
			end
			res ? Pass(:test,$(Expr(:quote, exprs)), nothing, nothing) :
		   	   Fail(:test,$(Expr(:quote, exprs)), fail_data, nothing)
		catch err
			Error(:test_error,$(Expr(:quote, exprs)), err, catch_backtrace())
		end end #quote #try
	else 		## logging is a path to file where the log is written to
		inner_ex = quote try
			log_file = open($logging,"a")
			writedlm(log_file,reshape([string(s) for s in $(var_data[:,1])],(1,$num_of_vars)),",")
			for n in 1:$maxtests
				$generate_values
				if $(esc(cond))
					num_good_args += 1
					res = res && $(esc(prop))
					writedlm(log_file,transpose(push!(convert(Vector{Any},$values),res)),",")
					if !res
						fail_data = $(esc(prop))
						break
					end
					if num_good_args >= $ntests
						break
					end
				end
			end
			close(log_file)
		        # no counter-example but insufficient amount of tests
			if res && num_good_args<$ntests
				nt = $ntests
 				error("Found only $num_good_args/$nt values satisfying given condition.")
			end
			res ? Pass(:test,$(Expr(:quote, exprs)), nothing, nothing) :
		   	   Fail(:test,$(Expr(:quote, exprs)), fail_data, nothing)
		catch err
			Error(:test_error,$(Expr(:quote, exprs)), err, catch_backtrace())
		end end #quote #try
	end
	inner_ex
#####
   catch err
	inner_ex = quote
		Error(:test_error,$(Expr(:quote, exprs)), $err, catch_backtrace())
	end
   end # try defining inner_ex
	return quote
		num_good_args = 0
		res = true
		fail_data = false
		result = $inner_ex
		record(get_testset(), result)
	end
end


text:
# export GetTheColor, printDict

# ======================================================================================
# EX: GetTheColor(border["color"])
# Get color depending on format and add opacity to it.
# ======================================================================================
function GetTheColor(shape, DOMColor)
                if isa(DOMColor, Array)
                    color = DOMColor
                # elseif isa(DOMColor, Float64)
                elseif isa(DOMColor, String)
                    if DOMColor[1:1] == "#"
                        if length(DOMColor) == 7
                            color = [ parse(Int, DOMColor[2:3], base = 16)*0.00390625,
                              parse(Int, DOMColor[4:5], base = 16)*0.00390625,
                              parse(Int, DOMColor[6:7], base = 16)*0.00390625 ]
                        elseif length(DOMColor) == 4
                            color = [ parse(Int, DOMColor[2:2]^2, base = 16)*0.00390625,
                              parse(Int, DOMColor[3:3]^2, base = 16)*0.00390625,
                              parse(Int, DOMColor[4:4]^2, base = 16)*0.00390625 ]
                        end
                    else
                        color = collect(color_names[ DOMColor ]).*0.00390625
                    end
                    # println("Color is $DOMColor")
                end

                if shape.flags[HasOpacity] && length(color) == 3
                    color = [color... , shape.opacity]
                    #push!(color, shape.opacity)
                    # println("opacity: ", typeof(color))
                    # println("opacity: ", color)
                end

    return color
end
# ======================================================================================
# Print out Element's DOM but not children
# Print out Element's DOM but not children
# ======================================================================================
function printDict(DOM)

    dict = copy(DOM)
    dict["nodes"] = "[...]"
    #contents = []

       keyList = sort(collect(keys(dict)))
       str, key, value = "H","",""
           for k in 1:length(keyList)
             key = keyList[k]
               if isa(dict[key], Dict)
                 value = "{ $(printDict(dict[key])) } "
               else
                 value = dict[keyList[k]]
               end
               if k != 1; key = ", $(key)"; end
            str =   "$(str)$(key):$(value)"
           end
    println(str)
    # return str
end

include("parseSml.jl")
include("ElementDefaults.jl")
include("DomTree.jl")
include("DomBuildDefault.jl")
include("DomShadow.jl")
include("DomToLayout.jl")


text:
@testset "Reducing ODE mod p" begin

    ode = @ODEmodel(
        x'(t) = 1//2 * x(t) - 5 * y(t),
        y'(t) = (3//5 * x(t) + 17 * y(t)) // (y(t) - 1),
        z(t) = y(t)
    )

    ode_red = reduce_ode_mod_p(ode, 17)

    x = str_to_var("x", ode_red.poly_ring)
    y = str_to_var("y", ode_red.poly_ring)
    @test ode_red.x_equations[x] == 9 * x + 12 * y
    @test ode_red.x_equations[y] == 4 * x // (y + 16)

end


text:
function init_gratis_start(context::PSCOPFContext, ref_schedule_type::DeciderType)

    if is_market(ref_schedule_type)
        reference_schedule = get_market_schedule(context)
    elseif is_tso(ref_schedule_type)
        reference_schedule = get_tso_schedule(context)
    else
        @sprintf("Invalid reference schedule type config : %s.", ref_schedule_type)
        throw( error(msg) )
    end
    gratis_starts = get_starts(reference_schedule, get_generators_initial_state(context))

    return gratis_starts
end

function get_starts(commitments::SortedDict{Tuple{String, Dates.DateTime}, GeneratorState}, initial_state::SortedDict{String, GeneratorState})
    result = Set{Tuple{String,Dates.DateTime}}()

    preceding_id, preceding_state = nothing, nothing
    for ((gen_id,ts), gen_state) in commitments
        if ( isnothing(preceding_id) || gen_id!=preceding_id )
            preceding_state = initial_state[gen_id]
        end

        if get_start_value(preceding_state, gen_state) > 1e-09
            push!(result, (gen_id,ts) )
        end

        preceding_id = gen_id
        preceding_state = gen_state
    end

    return result
end

function get_starts(schedule::Schedule, initial_state::SortedDict{String, GeneratorState})
    commitments = SortedDict{Tuple{String, Dates.DateTime}, GeneratorState}()

    for (gen_id, gen_schedule) in schedule.generator_schedules
        if isempty(gen_schedule.commitment)
            continue
        end
        for (ts, current_state) in gen_schedule.commitment
            if !is_definitive(current_state)
                break #if current state is not definitive the following starts cannot be definitive
            else
                commitments[gen_id, ts] = get_value(current_state)
            end
        end
    end

    return get_starts(commitments, initial_state)
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    @sampler(def)

Define a sampler type for some `scenariotype` compatible with StochasticPrograms using the syntax
```julia
@sampler samplername = begin
    ...internals...

    @sample scenariotype begin
        ...
        return scenario
    end
end
```
Any internal state required by the sampler, as well as any specialized constructor, are defined in the @sampler block as they would be in any Julia struct. Define the sample operation inside the [`@sample`](@ref) block and specify the `scenariotype` that the sampler returns. The defined sampler will be available on all Julia processes.

## Examples

The following defines a simple dummy sampler, with some internal weight value, for the scenario defined in [`@scenario`](@ref), and samples one scenario.

```jldoctest; setup = :(@scenario ExampleScenario = ξ::Float64), filter = r".*"
@sampler ExampleSampler = begin
    w::Float64

    ExampleSampler(w::AbstractFloat) = new(w)

    @sample ExampleScenario begin
        @parameters w
        return ExampleScenario(w*randn(), probability = rand())
    end
end
s = ExampleSampler(2.0)
s()

# output

ExampleScenario with probability 0.29
  ξ: 1.48

```

See also: [`@sample`](@ref), [`@scenario`](@ref)
"""
macro sampler(arg)
    @capture(prettify(arg), samplername_Symbol = samplerdef_) || error("Invalid syntax. Expected: sampler = begin samplerdef end")
    scenarioname = :undefined
    sampledefs = Vector{Expr}()
    def = postwalk(prettify(samplerdef)) do x
        if @capture(x, @parameters args__)
            code = Expr(:block)
            for param in args
                push!(code.args, :($param = sampler.$param))
            end
            return code
        elseif @capture(x, @sample sname_ sampledef_)
            scenarioname = sname
            push!(sampledefs, sampledef)
            return @q begin end
        else
            return x
        end
    end
    # Handle sample definition
    length(sampledefs) == 1 || error("Must provide exactly one @sample operation.")
    sampledef = sampledefs[1]
    # Define sampler type
    code = @q begin
        struct $samplername <: AbstractSampler{$scenarioname}
            $def
        end
        function (sampler::$samplername)()
            $sampledef
        end
    end
    code = prettify(code)
    code = @q begin
        if (@__MODULE__) == $(esc(:Main))
            @everywhere begin
                $code
            end
        else
            $(esc(code))
        end
    end
    return prettify(code)
end
"""
    @sample(def)

Define the sample operaton inside a @sampler block, using the syntax
```julia
@sample begin
    ...
    return sampled_scenario
end
```
The sampler object is referenced through the reserved keyword `sampler`, from which any internals can be accessed.
"""
macro sample(def) @warn "@sample should be used inside a @sampler block." end


text:
InputFolder = './Images/TwoPureColors/';
OutputFolder = './Results/Images/TwoPureColors/';

@iA = '*.tif';

@fxc_cggSplitTwoColorStains [iA] > [iNorm iH iE];
params.Io = 240;
params.beta = 0.15;
params.alpha = 1;
/endf

/show iNorm >;
/keep iNorm > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
include("euler/euler.jl")
using .Primes: prime_numbers

function compute(n::Int)::Int
    result = 1
    for prime ∈ prime_numbers(trunc(Int, √n))
        result *= prime
        if result * prime > n
            return result
        end
    end
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
struct ShaderSpecification
    source_file::String
    reuse_descriptors::Bool
    entry_point::Symbol
    stage::Vk.ShaderStageFlag
    language::ShaderLanguage
end

function ShaderSpecification(source_file, stage::Vk.ShaderStageFlag; reuse_descriptors = false, entry_point = :main)
    ShaderSpecification(source_file, reuse_descriptors, entry_point, stage, shader_language(source_file))
end

function ShaderSpecification(source_file, language::ShaderLanguage; reuse_descriptors = false, entry_point = :main)
    stage = @match language begin
        &GLSL || &HLSL => shader_stage(source_file, language)
        &SPIR_V => shader_stage(source_file, language) # will error, need to introspect into SPIR-V module
    end
    ShaderSpecification(source_file, reuse_descriptors, entry_point, stage, language)
end

function ShaderSpecification(source_file; reuse_descriptors = false, entry_point = :main)
    if shader_language(source_file) == SPIR_V
        ShaderSpecification(source_file, SPIR_V; reuse_descriptors, entry_point)
    else
        error("Language or stage must be supplied")
    end
end


text:
"""
    predict(branch,trunk,initial_condition,x_locations,t_values)

Predict solution ``u(t,x)`` at specified output locations using trained operator neural network `branch` and `trunk`.

"""
function predict(branch,trunk,initial_condition,x_locations,t_values)
    u = zeros(size(t_values,1),size(x_locations,1));
    bk = branch(initial_condition)';
    for i in 1:size(t_values,1)
        t_mat = repeat([t_values[i]],1,size(x_locations,1));
        t_x_mat = vcat(t_mat,x_locations');
        u[i,:] = bk*trunk(t_x_mat);
    end

    return u
end

"""
    loss_all(branch,trunk,initial_conditon,solution_location,target_value)

Compute the mean squared error (MSE) for a complete dataset.

"""
function loss_all(branch,trunk,initial_condition,solution_location,target_value)
    yhat = zeros(1,size(target_value,2));
    for i in 1:size(target_value,2)
        yhat[i] = branch(initial_condition[:,i])'*trunk(solution_location[:,i]);
    end
    return (1/size(target_value,2))*sum((yhat.-target_value).^2);
end

"""
    function build_dense_model(number_layers,neurons,activations)

Build a feedforward neural network (FFNN) consisting of `number_layers` of Flux dense layers for the specified number of `neurons` and `activations`.

# Examples
```julia-repl
julia> build_dense_model(2,[(128,128),(128,128)],[relu,relu])
Chain(Dense(128, 128, NNlib.relu), Dense(128, 128, NNlib.relu))
```

"""
function build_dense_model(number_layers,neurons,activations)
    layers = [Dense(neurons[i][1],neurons[i][2],activations[i]) for i in 1:number_layers];
    return Chain(layers...)|>f64
end

"""
    train_model(branch,trunk,n_epoch,train_data;learning_rate=0.00001,save_at=2500,starting_epoch=0)

Train the operator neural network using the mean squared error (MSE) and Adam optimization for `n_epochs` epochs.

"""
function train_model(branch,trunk,n_epoch,train_data,test_data,pde_function;learning_rate=1e-5,save_at=2500,starting_epoch=0)
    loss(x,y,z) = Flux.mse(branch(x)'*trunk(y),z)
    par = Flux.params(branch,trunk);
    opt = ADAM(learning_rate);
    @showprogress 1 "Training the model..." for i in 1:n_epoch
        Flux.train!(loss,par,train_data,opt);
        if i%save_at == 0
            save_model(branch,trunk,Int(starting_epoch+i),(pde_function*"_temp"))
            train_MSE = loss_all(branch,trunk,train_data.data[1],train_data.data[2],train_data.data[3]);
            test_MSE = loss_all(branch,trunk,test_data.data[1],test_data.data[2],test_data.data[3]);
            println("Train MSE $train_MSE")
            println("Test MSE $test_MSE")
        end
    end
    return branch,trunk
end

"""
    function exp_kernel_periodic(fnc,x_locations;length_scale=0.5)

Covariance kernel for radial basis function (GRF) and IC function for domain.

"""
function exp_kernel_periodic(fnc,x_locations,length_scale)
    out = zeros(size(x_locations,1),size(x_locations,1));
    out = [-(1/(2*length_scale^2))*norm(fnc(x_locations[i]) - fnc(x_locations[j]))^2 for i in 1:size(x_locations,1), j in 1:size(x_locations,1)]
    return exp.(out)
end

"""
    generate_periodic_functions(fnc,x_locations,number_functions,length_scale)

Generate a specified `number_functions` of random periodic vectors using the `exp_kernel_periodic` function and a multivariate distribution.

"""
function generate_periodic_functions(fnc,x_locations,number_functions,length_scale)
    sigma = exp_kernel_periodic(fnc,x_locations,length_scale); # Covariance
    mu = zeros(size(x_locations,1)); # Zero mean
    # Force the covariance matrix to be positive definite, by construction it is "approximately" Hermitian but Julia is strict
    mineig = eigmin(sigma);
    sigma -= mineig * I;
    d = MvNormal(mu,Symmetric(sigma)); # Multivariate distribution
    return rand(d,Int(number_functions))
end

"""
    generate_sinusoidal_functions_2_parameter(x_locations,number_functions)

Generate a specified `number_functions` of random periodic vectors for the distribution \$α \\sin(x)+β\$ for ``α ∈ [-1,1]`` and ``β ∈ [-1,1]``.

"""
function generate_sinusoidal_functions_2_parameter(x_locations,number_functions)
    values = rand(-1.0:(2.0/(number_functions*10)):1.0,(number_functions,2)); # alpha, beta parameters
    initial_base = vcat(values,[0.0 0.0]);
    if size(unique(initial_base,dims=2),1) != size(initial_base,1)
        error("Duplicate multiplier sets or generating function included in data!")
    end
    initial_conditions = zeros(size(x_locations,1),number_functions);
    for i in 1:number_functions
        initial_conditions[:,i] = values[i,1]*sin.(x_locations).+values[i,2]; # αsin(x)+β
    end
    return initial_conditions
end

"""
    solution_extraction(x_locations,t_values,solution,initial_condition,number_solution_points)

Extract the specified `number_solution_points` randomly from the ``u(t,x)`` solution space.

"""
function solution_extraction(x_locations,t_values,solution,initial_condition,number_solution_points)
    if number_solution_points >= size(solution,1)*size(solution,2)
        error("Invalid number of test points for the given dataset size!")
    else
        # Create a grid of the t,x data in the form of an array of tuples for accessing below
        t_x_grid = Array{Tuple{Float64,Float64}}(undef,(size(solution,1), size(solution,2)));
        for i in 1:size(solution,1)
            for j in 1:size(solution,2)
                t_x_grid[i,j] = (t_values[i],x_locations[j]);
            end
        end
        lin_ind = reduce(vcat,size(solution[:]));
        shuffled_indices = randperm(lin_ind) # Randomly shuffle the linear indices corresponding to the t,x data
        indices = shuffled_indices[1:number_solution_points] # From shuffled data, extract the number of test points for training
        initial = repeat(reshape(initial_condition,(length(initial_condition),1)),1,number_solution_points);
        sol_location = zeros(2,number_solution_points);
        for i in 1:2
            for j in 1:number_solution_points
                sol_location[i,j] = t_x_grid[indices[j]][i]
            end
        end
        sol = reshape(solution[indices],(1,number_solution_points));
        return initial, sol_location, sol
    end
end

"""
    generate_periodic_train_test_initial_conditions(t_span,number_sensors,number_test_functions,number_train_functions,number_solution_points,pde_function_handle;L1=0,L2=2*pi,length_scale=0.5,batch=number_solution_points,dt=1e-3,nu_val=0.1,domain="periodic",fnc=(x)->sin(x/2)^2)

Generate the training and testing data for a specified `pde_function_handle` for periodic boundary conditions using a Fourier spectral method. Defaults to IC \$f(\\sin^2(x/2))\$ and \$x ∈ [0,1]\$.

"""
function generate_periodic_train_test_initial_conditions(number_sensors,number_train_functions,number_test_functions;L1=0,L2=2*pi,length_scale=0.5,domain="periodic",fnc=(x)->sin(x/2)^2)

    if domain == "periodic"
        x_full = range(L1,stop = L2,length = number_sensors+1); # Full domain for periodic number of sensors
        random_ics = generate_periodic_functions(fnc,x_full,Int(number_train_functions + number_test_functions),length_scale);
    elseif domain == "full"
        x_full = range(L1,stop = L2,length = number_sensors); # Full domain for periodic number of sensors
        random_ics = generate_periodic_functions(fnc,x_full,Int(number_train_functions + number_test_functions),length_scale);
    end

    train_ic = random_ics[:,1:number_train_functions];
    test_ic = random_ics[:,number_train_functions+1:end];

    return train_ic, test_ic
end

"""
    generate_periodic_train_test(t_span,number_sensors,number_test_functions,number_train_functions,number_solution_points,pde_function_handle;L1=0,L2=2*pi,length_scale=0.5,batch=number_solution_points,dt=1e-3,nu_val=0.1,domain="periodic",fnc=(x)->sin(x/2)^2)

Generate the training and testing data for a specified `pde_function_handle` for periodic boundary conditions using a Fourier spectral method. Defaults to IC \$f(\\sin^2(x/2))\$ and \$x ∈ [0,1]\$.

"""
function generate_periodic_train_test(t_span,number_sensors,number_train_functions,number_test_functions,number_solution_points,pde_function_handle;L1=0,L2=2*pi,length_scale=0.5,batch=number_solution_points,dt_size=1e-3,nu_val=0.1,domain="periodic",fnc=(x)->sin(x/2)^2)

    if domain == "periodic"
        x_full = range(L1,stop = L2,length = number_sensors+1); # Full domain for periodic number of sensors
        random_ics = generate_periodic_functions(fnc,x_full,Int(number_train_functions + number_test_functions),length_scale)
        dL = abs(L2-L1);
        # Set up x domain and wave vector for spectral solution
        x = trapezoid(number_sensors,L1,L2)[1];
        k = reduce(vcat,(2*π/dL)*[0:number_sensors/2-1 -number_sensors/2:-1]);
    elseif domain == "full"
        x_full = range(L1,stop = L2,length = number_sensors); # Full domain for periodic number of sensors
        random_ics = generate_periodic_functions(fnc,x_full,Int(number_train_functions + number_test_functions),length_scale)
        dL = abs(L2-L1);
        # Set up x domain and wave vector for spectral solution
        x = trapezoid((number_sensors-1),L1,L2)[1];
        k = reduce(vcat,(2*π/dL)*[0:Int(number_sensors-1)/2-1 -Int(number_sensors-1)/2:-1]);
    end

    # Generate the dataset using spectral method
    t_length = t_span[2]/dt_size + 1;
    t = range(t_span[1],stop = t_span[2], length = Int(t_length));
    train_ic = zeros(number_sensors,number_solution_points,number_train_functions);
    train_loc = zeros(2,number_solution_points,number_train_functions);
    train_target = zeros(1,number_solution_points,number_train_functions);
    test_ic = zeros(number_sensors,number_solution_points,number_test_functions);
    test_loc = zeros(2,number_solution_points,number_test_functions);
    test_target = zeros(1,number_solution_points,number_test_functions);

    if domain == "periodic"
        @showprogress 1 "Building training dataset..." for i in 1:number_train_functions
            interp_train_sample = random_ics[1:end-1,i];
            u_train = generate_fourier_solution(L1,L2,t_span,number_sensors,interp_train_sample,pde_function_handle,dt=dt_size)[1];
            train_ic[:,:,i], train_loc[:,:,i], train_target[:,:,i] = solution_extraction(x,t,u_train,interp_train_sample,number_solution_points);
        end
        @showprogress 1 "Building testing dataset..." for i in 1:number_test_functions
            interp_test_sample = random_ics[1:end-1,Int(number_train_functions+i)];
            u_test = generate_fourier_solution(L1,L2,t_span,number_sensors,interp_test_sample,pde_function_handle,dt=dt_size)[1];
            test_ic[:,:,i], test_loc[:,:,i], test_target[:,:,i] = solution_extraction(x,t,u_test,interp_test_sample,number_solution_points);
        end
    elseif domain == "full"
        @showprogress 1 "Building training dataset..." for i in 1:number_train_functions
            interp_train_sample_fourier = random_ics[1:end-1,i];
            u_train_fourier = generate_fourier_solution(L1,L2,t_span,Int(number_sensors-1),interp_train_sample_fourier,pde_function_handle,dt=dt_size)[1];
            interp_train_sample = random_ics[:,i];
            u_train = periodic_fill_solution(u_train_fourier);
            train_ic[:,:,i], train_loc[:,:,i], train_target[:,:,i] = solution_extraction(x_full,t,u_train,interp_train_sample,number_solution_points);
        end
        @showprogress 1 "Building testing dataset..." for i in 1:number_test_functions
            interp_test_sample_fourier = random_ics[1:end-1,Int(number_train_functions+i)];
            u_test_fourier = generate_fourier_solution(L1,L2,t_span,Int(number_sensors-1),interp_test_sample_fourier,pde_function_handle,dt=dt_size)[1];
            interp_test_sample = random_ics[:,Int(number_train_functions+i)];
            u_test = periodic_fill_solution(u_test_fourier);
            test_ic[:,:,i], test_loc[:,:,i], test_target[:,:,i] = solution_extraction(x_full,t,u_test,interp_test_sample,number_solution_points);
        end
    end

    # Combine data sets from each function
    opnn_train_ic = reshape(hcat(train_ic...),(number_sensors,Int(number_solution_points*number_train_functions)));
    opnn_train_loc = reshape(hcat(train_loc...),(2,Int(number_solution_points*number_train_functions)));
    opnn_train_target = reshape(hcat(train_target...),(1,Int(number_solution_points*number_train_functions)));
    opnn_test_ic = reshape(hcat(test_ic...),(number_sensors,Int(number_solution_points*number_test_functions)));
    opnn_test_loc = reshape(hcat(test_loc...),(2,Int(number_solution_points*number_test_functions)));
    opnn_test_target = reshape(hcat(test_target...),(1,Int(number_solution_points*number_test_functions)));

    train_data = DataLoader(opnn_train_ic, opnn_train_loc, opnn_train_target, batchsize=batch);
    test_data = DataLoader(opnn_test_ic, opnn_test_loc, opnn_test_target, batchsize=batch);

    return train_data, test_data
end

"""
    generate_periodic_train_test_initial_condition_load(train_ics,test_ics,t_span,number_sensors,number_test_functions,number_train_functions,number_solution_points,pde_function_handle;L1=0,L2=2*pi,batch=number_solution_points,dt=1e-3,nu_val=0.1,domain="periodic")

Generate the training and testing data for a specified `pde_function_handle` for periodic boundary conditions using a Fourier spectral method. Loads previously generated initial conditions, see [`generate_periodic_train_test`](@ref) for generating all data in one go.

"""
function generate_periodic_train_test_initial_condition_load(train_ics,test_ics,t_span,number_sensors,number_train_functions,number_test_functions,number_solution_points,pde_function_handle;L1=0,L2=2*pi,batch=number_solution_points,dt_size=1e-3,nu_val=0.1,domain="periodic")

    if domain == "periodic"
        dL = abs(L2-L1);
        # Set up x domain and wave vector for spectral solution
        x = trapezoid(number_sensors,L1,L2)[1];
        k = reduce(vcat,(2*π/dL)*[0:number_sensors/2-1 -number_sensors/2:-1]);
    elseif domain == "full"
        dL = abs(L2-L1);
        # Set up x domain and wave vector for spectral solution
        x = trapezoid((number_sensors-1),L1,L2)[1];
        k = reduce(vcat,(2*π/dL)*[0:Int(number_sensors-1)/2-1 -Int(number_sensors-1)/2:-1]);
    end

    # Generate the dataset using spectral method
    t_length = t_span[2]/dt_size + 1;
    t = range(t_span[1],stop = t_span[2], length = Int(t_length));
    train_ic = zeros(number_sensors,number_solution_points,number_train_functions);
    train_loc = zeros(2,number_solution_points,number_train_functions);
    train_target = zeros(1,number_solution_points,number_train_functions);
    test_ic = zeros(number_sensors,number_solution_points,number_test_functions);
    test_loc = zeros(2,number_solution_points,number_test_functions);
    test_target = zeros(1,number_solution_points,number_test_functions);

    if domain == "periodic"
        @showprogress 1 "Building training dataset..." for i in 1:number_train_functions
            interp_train_sample = train_ics[1:end-1,i];
            u_train = generate_fourier_solution(L1,L2,t_span,number_sensors,interp_train_sample,pde_function_handle,dt=dt_size)[1];
            train_ic[:,:,i], train_loc[:,:,i], train_target[:,:,i] = solution_extraction(x,t,u_train,interp_train_sample,number_solution_points);
        end
        @showprogress 1 "Building testing dataset..." for i in 1:number_test_functions
            interp_test_sample = test_ics[1:end-1,i];
            u_test = generate_fourier_solution(L1,L2,t_span,number_sensors,interp_test_sample,pde_function_handle,dt=dt_size)[1];
            test_ic[:,:,i], test_loc[:,:,i], test_target[:,:,i] = solution_extraction(x,t,u_test,interp_test_sample,number_solution_points);
        end
    elseif domain == "full"
        @showprogress 1 "Building training dataset..." for i in 1:number_train_functions
            interp_train_sample_fourier = train_ics[1:end-1,i];
            u_train_fourier = generate_fourier_solution(L1,L2,t_span,Int(number_sensors-1),interp_train_sample_fourier,pde_function_handle,dt=dt_size)[1];
            interp_train_sample = random_ics[:,i];
            u_train = periodic_fill_solution(u_train_fourier);
            train_ic[:,:,i], train_loc[:,:,i], train_target[:,:,i] = solution_extraction(x_full,t,u_train,interp_train_sample,number_solution_points);
        end
        @showprogress 1 "Building testing dataset..." for i in 1:number_test_functions
            interp_test_sample_fourier = test_ics[1:end-1,i];
            u_test_fourier = generate_fourier_solution(L1,L2,t_span,Int(number_sensors-1),interp_test_sample_fourier,pde_function_handle,dt=dt_size)[1];
            interp_test_sample = random_ics[:,Int(number_train_functions+i)];
            u_test = periodic_fill_solution(u_test_fourier);
            test_ic[:,:,i], test_loc[:,:,i], test_target[:,:,i] = solution_extraction(x_full,t,u_test,interp_test_sample,number_solution_points);
        end
    end

    # Combine data sets from each function
    opnn_train_ic = reshape(hcat(train_ic...),(number_sensors,Int(number_solution_points*number_train_functions)));
    opnn_train_loc = reshape(hcat(train_loc...),(2,Int(number_solution_points*number_train_functions)));
    opnn_train_target = reshape(hcat(train_target...),(1,Int(number_solution_points*number_train_functions)));
    opnn_test_ic = reshape(hcat(test_ic...),(number_sensors,Int(number_solution_points*number_test_functions)));
    opnn_test_loc = reshape(hcat(test_loc...),(2,Int(number_solution_points*number_test_functions)));
    opnn_test_target = reshape(hcat(test_target...),(1,Int(number_solution_points*number_test_functions)));

    train_data = DataLoader(opnn_train_ic, opnn_train_loc, opnn_train_target, batchsize=batch);
    test_data = DataLoader(opnn_test_ic, opnn_test_loc, opnn_test_target, batchsize=batch);

    return train_data, test_data
end

"""
    generate_periodic_train_test_muscl(L1,L2,t_span,number_sensors,number_test_functions,number_train_functions,number_solution_points,pde_function_handle;L1=0,L2=2*pi,length_scale=0.5,batch=number_solution_points,dt_size=1e-4,upwind_solution_points=4096,fnc=(x)->sin(x/2)^2)

Generate the training and testing data for a specified `pde_function_handle` for periodic boundary conditions using a MUSCL method. Defaults to IC \$f(\\sin^2(x/2))\$ and \$x ∈ [0,1]\$.

"""
function generate_periodic_train_test_muscl(t_span,number_sensors,number_train_functions,number_test_functions,number_solution_points,pde_function_handle;L1=0,L2=2*pi,length_scale=0.5,batch=number_solution_points,dt_size=1e-4,upwind_solution_points=4096,fnc=(x)->sin(x/2)^2)

    x_full = range(L1,stop = L2,length = upwind_solution_points+1); # Full domain
    random_ics = generate_periodic_functions(fnc,x_full,Int(number_train_functions + number_test_functions),length_scale)
    dL = abs(L2-L1);

    # Set up x domains
    x = trapezoid(number_sensors,L1,L2)[1];
    x_full = trapezoid(upwind_solution_points,L1,L2)[1];

    # Generate the dataset
    t_length = t_span[2]/dt_size + 1;
    t = range(t_span[1],stop = t_span[2], length = Int(t_length));
    train_ic = zeros(Int(number_sensors),number_solution_points,number_train_functions);
    train_loc = zeros(2,number_solution_points,number_train_functions);
    train_target = zeros(1,number_solution_points,number_train_functions);
    test_ic = zeros(Int(number_sensors),number_solution_points,number_test_functions);
    test_loc = zeros(2,number_solution_points,number_test_functions);
    test_target = zeros(1,number_solution_points,number_test_functions);

    @showprogress 1 "Building training dataset..." for i in 1:number_train_functions

        interp_train_sample_full = random_ics[1:end-1,i];
        u_train_full = generate_muscl_minmod_solution(L1,L2,t_span[end],upwind_solution_points,interp_train_sample_full,muscl_minmod_RHS!;dt=dt_size)
        interp_train_sample = transpose(solution_spatial_sampling(x,x_full,transpose(interp_train_sample_full)));
        u_train = solution_spatial_sampling(x,x_full,u_train_full);
        train_ic[:,:,i], train_loc[:,:,i], train_target[:,:,i] = solution_extraction(x,t,u_train,interp_train_sample,number_solution_points);
    end

    @showprogress 1 "Building testing dataset..." for i in 1:number_test_functions

        interp_test_sample_full = random_ics[1:end-1,Int(number_train_functions+i)];
        u_test_full = generate_muscl_minmod_solution(L1,L2,t_span[end],upwind_solution_points,interp_test_sample_full,muscl_minmod_RHS!;dt=dt_size);
        interp_test_sample = transpose(solution_spatial_sampling(x,x_full,transpose(interp_test_sample_full)));
        u_test = solution_spatial_sampling(x,x_full,u_test_full);
        test_ic[:,:,i], test_loc[:,:,i], test_target[:,:,i] = solution_extraction(x,t,u_test,interp_test_sample,number_solution_points);
    end

    # Combine data sets from each function
    opnn_train_ic = reshape(hcat(train_ic...),(Int(number_sensors),Int(number_solution_points*number_train_functions)));
    opnn_train_loc = reshape(hcat(train_loc...),(2,Int(number_solution_points*number_train_functions)));
    opnn_train_target = reshape(hcat(train_target...),(1,Int(number_solution_points*number_train_functions)));
    opnn_test_ic = reshape(hcat(test_ic...),(Int(number_sensors),Int(number_solution_points*number_test_functions)));
    opnn_test_loc = reshape(hcat(test_loc...),(2,Int(number_solution_points*number_test_functions)));
    opnn_test_target = reshape(hcat(test_target...),(1,Int(number_solution_points*number_test_functions)));

    train_data = DataLoader(opnn_train_ic, opnn_train_loc, opnn_train_target, batchsize=batch);#,shuffle = true);
    test_data = DataLoader(opnn_test_ic, opnn_test_loc, opnn_test_target, batchsize=batch);
    return train_data, test_data
end

"""
    save_model(branch,trunk,n_epoch,loss_all_train,loss_all_test,pde_function)

Save the trained `branch` and `trunk` neural networks and the training and testing loss history.

"""
function save_model(branch,trunk,n_epoch,pde_function)
    @save @sprintf("branch_epochs_%i_%s.bson",n_epoch,pde_function) branch
    @save @sprintf("trunk_epochs_%i_%s.bson",n_epoch,pde_function) trunk
end

"""
    load_model(n_epoch,pde_function)

Load the trained `branch` and `trunk` neural networks.

"""
function load_model(n_epoch,pde_function)
    @load @sprintf("branch_epochs_%i_%s.bson",n_epoch,pde_function) branch
    @load @sprintf("trunk_epochs_%i_%s.bson",n_epoch,pde_function) trunk
    return branch, trunk
end

"""
    load_branch(n_epoch,pde_function)

Load the trained `branch` neural networks.

"""
function load_branch(n_epoch,pde_function)
    @load @sprintf("branch_epochs_%i_%s.bson",n_epoch,pde_function) branch
    return branch
end

"""
    load_trunk(n_epoch,pde_function)

Load the trained `trunk` neural networks.

"""
function load_trunk(n_epoch,pde_function)
    @load @sprintf("trunk_epochs_%i_%s.bson",n_epoch,pde_function) trunk
    return trunk
end

"""
    save_data(train_data,test_data,number_train_functions,number_test_functions,number_solution_points,pde_function)

Save the `train_data` and `test_data`.

"""
function save_data(train_data,test_data,number_train_functions,number_test_functions,number_solution_points,pde_function)
    train_ic = train_data.data[1];
    train_loc = train_data.data[2];
    train_sol = train_data.data[3];
    test_ic = test_data.data[1];
    test_loc = test_data.data[2];
    test_sol = test_data.data[3];
    @save @sprintf("train_ic_data_%i_%s.bson",number_train_functions,pde_function) train_ic
    @save @sprintf("test_ic_data_%i_%s.bson",number_test_functions,pde_function) test_ic
    @save @sprintf("train_loc_data_%i_%s.bson",number_train_functions,pde_function) train_loc
    @save @sprintf("test_loc_data_%i_%s.bson",number_test_functions,pde_function) test_loc
    @save @sprintf("train_target_data_%i_%s.bson",number_train_functions,pde_function) train_sol
    @save @sprintf("test_target_data_%i_%s.bson",number_test_functions,pde_function) test_sol
end

"""
    load_data(n_epoch,number_train_functions,number_test_functions,pde_function)

Load the trained `branch` and `trunk` neural networks along with the `train_data` and `test_data`.

"""
function load_data(n_epoch,number_train_functions,number_test_functions,pde_function)
    @load @sprintf("branch_epochs_%i_%s.bson",n_epoch,pde_function) branch
    @load @sprintf("trunk_epochs_%i_%s.bson",n_epoch,pde_function) trunk
    @load @sprintf("train_ic_data_%i_%s.bson",number_train_functions,pde_function) train_ic
    @load @sprintf("train_loc_data_%i_%s.bson",number_train_functions,pde_function) train_loc
    @load @sprintf("train_target_data_%i_%s.bson",number_train_functions,pde_function) train_sol
    @load @sprintf("test_ic_data_%i_%s.bson",number_test_functions,pde_function) test_ic
    @load @sprintf("test_loc_data_%i_%s.bson",number_test_functions,pde_function) test_loc
    @load @sprintf("test_target_data_%i_%s.bson",number_test_functions,pde_function) test_sol
    return branch, trunk, train_ic, train_loc, train_sol, test_ic, test_loc, test_sol#
end

"""
    load_data_train_test(number_train_functions,number_test_functions,pde_function)

Load the `train_data` and `test_data`.

"""
function load_data_train_test(number_train_functions,number_test_functions,pde_function)
    @load @sprintf("train_ic_data_%i_%s.bson",number_train_functions,pde_function) train_ic
    @load @sprintf("train_loc_data_%i_%s.bson",number_train_functions,pde_function) train_loc
    @load @sprintf("train_target_data_%i_%s.bson",number_train_functions,pde_function) train_sol
    @load @sprintf("test_ic_data_%i_%s.bson",number_test_functions,pde_function) test_ic
    @load @sprintf("test_loc_data_%i_%s.bson",number_test_functions,pde_function) test_loc
    @load @sprintf("test_target_data_%i_%s.bson",number_test_functions,pde_function) test_sol
    return train_ic, train_loc, train_sol, test_ic, test_loc, test_sol
end

"""
    load_data_initial_conditions(number_train_functions,number_test_functions,pde_function)

Load the initial conditions from the `train_data` and `test_data`.

"""
function load_data_initial_conditions(number_train_functions,number_test_functions,pde_function)
    @load @sprintf("train_ic_data_%i_%s.bson",number_train_functions,pde_function) train_ic
    @load @sprintf("test_ic_data_%i_%s.bson",number_test_functions,pde_function) test_ic
    return train_ic, test_ic
end

"""
    save_data_initial_conditions(number_train_functions,number_test_functions)

Saves the initial conditions from the `train_ic` and `test_ic`.

"""
function save_data_initial_conditions(train_ic,test_ic,number_train_functions,number_test_functions,pde_function)
    @save @sprintf("train_ic_data_%i_%s.bson",number_train_functions,pde_function) train_ic
    @save @sprintf("test_ic_data_%i_%s.bson",number_test_functions,pde_function) test_ic
end


text:
function prob(n)
  n = BigInt(n)
  two = BigInt(2)
  x = div(n,2)
  a = factorial(n)
  b = factorial(x)
  b *= b
  c = two^n
  return a // b // c
end

for n in [3030, 1000, 500, 200, 100, 50, 10]
  p = float64(round(prob(n)*100, 2))
  println(n, " : ", p, " %")
end

max_stars_repo_name: jgoizueta/binomial-techniques

max_issues_repo_name: jgoizueta/binomial-techniques

max_forks_repo_name: jgoizueta/binomial-techniques


text:
function sample_full_factorial(factors::Array{T, 1}) where T <: Any
    return Array{Any, 1}([rand(i) for i in factors])
end

function check_repeated_row(subset::Array{T, 2}, row::Array{T, 1}) where T <: Any
    for j = 1:size(subset, 1)
        if subset[j, :] == row
            return true
        end
    end

    return false
end

function full_factorial_subset(factors::OrderedDict, experiments::Int)
    subset = fill!(Array{Any, 2}(experiments, length(keys(factors))), 0.0)

    for i = 1:experiments
        sample_row = sample_full_factorial(collect(values(factors)))

        while check_repeated_row(subset, sample_row)
            sample_row = sample_full_factorial(collect(values(factors)))
        end

        subset[i, :] = sample_row
    end

    subset = DataFrame(subset)
    rename!(subset, zip(names(subset), collect(keys(factors))))

    for factor in keys(factors)
        if eltype(subset[factor]) != eltype(factors[factor])
            subset[factor] = convert(Array{eltype(factors[factor]), 1}, subset[factor])
        end
    end

    return subset
end

function enforce_bounds(factors::OrderedDict,
                        sample_range::UnitRange{Int},
                        designs::Int,
                        check_bounds::Bool)
    println("> Factors: ", factors)

    full_factorial_size    = prod(length, values(factors))
    full_factorial_subsets = 2.0 ^ full_factorial_size

    println("> Full Factorial Size: ", full_factorial_size)
    println("> Total Subsets: ", full_factorial_subsets)
    println("> Range of Design Sizes: ", sample_range)
    println("> Number of Designs to Sample: ", designs)

    if check_bounds
        if sample_range.start == sample_range.stop
            restricted_subsets = factorial(float(full_factorial_size)) /
                                 (factorial(float(full_factorial_size - sample_range.start)) *
                                 factorial(float(sample_range.start)))
            println("> Total Subsets for Fixed Size ",
                    sample_range.start, ": ",
                    restricted_subsets)

            if designs > restricted_subsets
                println("> Requested too many designs, using ",
                        restricted_subsets, " instead")
                designs = floor(Int, restricted_subsets)
            end
        elseif designs > full_factorial_subsets
            println("> Requested too many designs, using ",
                    full_factorial_subsets, " instead")
            designs = floor(Int, full_factorial_subsets)
        end

        if sample_range.stop > full_factorial_size
            println("> Requested too many maximum experiments, using ",
                    full_factorial_size, " instead")
            sample_range = sample_range.start:floor(Int, full_factorial_size)
        end
    else
        println("> WARNING: Skipping bounds check!")
    end

    return designs, sample_range
end

function evaluate_all_metrics(factors::OrderedDict,
                              formula::Formula,
                              sample_range::UnitRange{Int},
                              designs::Int,
                              scale::Function)
    evaluation = DataFrame(Length  = [],
                           D       = [],
                           DELB    = [],
                           DELB_ad = [],
                           A       = [],
                           V       = [],
                           G       = [],
                           CN      = [],
                           GE      = [],
                           log2CN  = [],
                           log10D  = [])

    for i in 1:designs
        samples   = rand(sample_range)
        subset    = full_factorial_subset(factors, samples)
        candidate = generate_model_matrix(formula, subset, factors,
                                          scale = scale)

        d_opt = d_optimality(candidate)
        c_n   = condition_number(candidate)

        push!(evaluation, [size(candidate, 1),
                           d_opt,
                           d_efficiency_lower_bound(candidate),
                           d_efficiency_lower_bound_algdesign(candidate),
                           a_optimality(candidate),
                           v_optimality(candidate),
                           g_optimality(candidate),
                           c_n,
                           g_efficiency(candidate),
                           log(2, abs(c_n)),
                           log(10, abs(d_opt))])
    end

    return evaluation
end

function evaluate_metrics(factors::OrderedDict,
                          formula::Formula,
                          sample_range::UnitRange{Int},
                          designs::Int,
                          scale::Function)
    evaluation = DataFrame(Length  = [],
                           D       = [],
                           DELB    = [])

    for i in 1:designs
        samples   = rand(sample_range)
        subset    = full_factorial_subset(factors, samples)
        candidate = generate_model_matrix(formula, subset, factors,
                                          scale = scale)

        push!(evaluation, [size(candidate, 1),
                           d_optimality(candidate),
                           d_efficiency_lower_bound(candidate)])
    end

    return evaluation
end

function generate_designs(factors::OrderedDict,
                          formula::Formula,
                          sample_range::UnitRange{Int},
                          designs::Int;
                          check_bounds::Bool = true,
                          scale::Function = scale_boxdraper_encoding!,
                          compute_all_metrics::Bool = false)
    designs, sample_range = enforce_bounds(factors, sample_range,
                                           designs, check_bounds)

    if compute_all_metrics
        return evaluate_all_metrics(factors::OrderedDict,
                                    formula::Formula,
                                    sample_range::UnitRange{Int},
                                    designs::Int,
                                    scale::Function)
    else
        return evaluate_metrics(factors::OrderedDict,
                                formula::Formula,
                                sample_range::UnitRange{Int},
                                designs::Int,
                                scale::Function)
    end
end

function sample_subset(factors::OrderedDict,
                       sample_range::UnitRange{Int},
                       designs::Int;
                       check_bounds::Bool = true,
                       scale::Function = scale_boxdraper_encoding!)
    formula = build_linear_formula(collect(keys(factors)))

    # Any formula could potentially be used here:
    #   formula = @formula(y ~ x1 * x2 + x3)

    run_time = @elapsed sampling_subset = generate_designs(factors,
                                                           formula,
                                                           sample_range,
                                                           designs,
                                                           check_bounds = check_bounds,
                                                           scale = scale)
    println("> Elapsed Time: ", run_time, " seconds")

    sort!(sampling_subset, :D, rev = true)

    return sampling_subset
end

function sample_subsets(factors::Array{OrderedDict{Symbol, Any}, 1},
                        ranges::Array{UnitRange{Int}, 1},
                        designs::Int;
                        check_bounds::Bool = true,
                        scale::Function = scale_boxdraper_encoding!,
                        compute_all_metrics::Bool = false)
    sampled_subsets = []

    for subset = 1:length(ranges)
        sampled_subset = sample_subset(factors[subset],
                                       ranges[subset],
                                       designs,
                                       check_bounds = check_bounds,
                                       scale = scale)

        push!(sampled_subsets, sampled_subset)
    end

    return sampled_subsets
end


text:
using InducingPoints, KernelFunctions
using BenchmarkTools
using Plots
plotlyjs()

nv = [2, 4, 8, 16]
nnv = length(nv)
ker = SqExponentialKernel()

sizeug = zeros(nnv)
sizeexpgrid = zeros(nnv)

kermtimeug = zeros(nnv)
kermtimegrid = zeros(nnv)
kermmemug = zeros(nnv)
kermmemgrid = zeros(nnv)

sumtimeug = zeros(nnv)
sumtimegrid = zeros(nnv)
summemug = zeros(nnv)
summemgrid = zeros(nnv)

for (i, n) in enumerate(nv)
    iter = LinRange(0.0, 10.0, n)
    proditer = Iterators.product(iter, iter)
    ug = UniformGrid(proditer)
    expgrid = ug[:]
    println(length(expgrid))

    sizeug[i] = sizeof(ug)
    sizeexpgrid[i] = sizeof(expgrid)

    r1 = @benchmark kernelmatrix($ker, $ug)
    kermtimeug[i] = median(r1.times)
    kermmemug[i] = r1.memory
    r2 = @benchmark kernelmatrix($ker, $expgrid)
    kermtimegrid[i] = median(r2.times)
    kermmemgrid[i] = r2.memory

    rs1 = @benchmark sum.($ug)
    sumtimeug[i] = median(rs1.times)
    summemug[i] = rs1.memory
    rs2 = @benchmark sum.($expgrid)
    sumtimegrid[i] = median(rs2.times)
    summemgrid[i] = rs2.memory
end

elbs = ["" ""]
p1 = plot()
plot!(
    p1,
    nv,
    [sizeug, sizeexpgrid] ./ 1e3;
    label=["UniformGrid" "Explicit Grid"],
    xticks=nv,
    linewidth=3,
    legend=:topleft,
    xlabel="Nr Grid Points",
    ylabel="Object Size [kB]",
)

p2 = plot()
plot!(
    nv,
    [kermtimeug, kermtimegrid] ./ 1e3;
    yaxis=:log,
    xaxis=:log,
    label=elbs,
    title="Kernelmatrix",
    ylabel="Compute Time [μs]",
    left_margin=5.0Plots.Measures.mm,
    xticks=(nv, nv),
    linewidth=3,
)

p3 = plot()
plot!(
    nv,
    [sumtimeug, sumtimegrid] ./ 1e3;
    yaxis=:log,
    xaxis=:log,
    label=elbs,
    title="Broadcast",
    xticks=(nv, nv),
    linewidth=3,
    left_margin=-2.0Plots.Measures.mm,
)

p4 = plot()
plot!(
    p4,
    nv,
    [kermmemug, kermmemgrid] ./ 1e3;
    yaxis=:log,
    xaxis=:log,
    label=elbs,
    xlabel="Nr Grid Points",
    ylabel="Memory [kB]",
    left_margin=5.0Plots.Measures.mm,
    xticks=(nv, nv),
    linewidth=3,
)

p5 = plot()
plot!(
    p5,
    nv,
    [summemug, summemgrid] ./ 1e3;
    yaxis=:log,
    xaxis=:log,
    label=elbs,
    xlabel="Nr of Grid Points",
    left_margin=-2.0Plots.Measures.mm,
    xticks=(nv, nv),
    linewidth=3,
)

pt = plot(p2, p3)
pm = plot(p4, p5)
pR = plot(pt, pm; layout=(2, 1))

p = plot(
    p1, pR; size=(700, 400), margin=3Plots.Measures.mm, layout=grid(1, 2; widths=[0.3, 0.7])
)

savefig(p, "../src/assets/UniformGrid_bench.svg")


text:
using SparseArrays: SparseMatrixCSC, sparse
using LinearAlgebra: I
using OSQP: Ccsc, ManagedCcsc
@testset "sparse matrix interface roundtrip" begin
    jl = sparse(Matrix{Bool}(LinearAlgebra.I, 5, 5))
    mc = ManagedCcsc(jl)
    c = Ccsc(mc)
    jl2 = convert(SparseMatrixCSC, c)
    @test jl == jl2
end