text:
# abstract FieldState
# type NotInField          <: FieldState end
# type InQuotedField       <: FieldState end
# type InUnQuotedField     <: FieldState end
# type InDoubleQuotedField <: FieldState end

# not_in_field           = NotInField()
# in_quoted_field        = InQuotedField()
# in_un_quoted_field     = InUnQuotedField()
# in_double_quoted_field = InDoubleQuotedField()

baremodule FieldState
    not_in_field           = 1
    in_quoted_field        = 2
    in_un_quoted_field     = 3
    in_double_quoted_field = 4
end

function get_string_type(strtype::String)
	if lowercase(strtype)=="ascii"
		return ASCIIString
	elseif lowercase(strtype)=="utf8"
		return UTF8String
	end
	throw(@sprintf("Invalid String Type %s", strtype))
end

function _iterxsv(io::IO; delimiter=','::Char, quotechar='"'::Char, strtype="ascii", close_io=false::Bool)
    string_type = get_string_type(strtype)
    state = FieldState.not_in_field
    delimiter_n = uint8(delimiter)
    quotechar_n = uint8(quotechar)
    r_newline_n = 0x0d
    n_newline_n = 0x0a
    row     = Array(string_type, 0)
    field   = Array(Uint8, 1024)
    n_field = 0

    while !eof(io)
        c = read(io, Uint8)
        if state==FieldState.not_in_field
            if c==quotechar_n
                state = FieldState.in_quoted_field
            elseif c==r_newline_n || c==n_newline_n
                # state remains the same
            elseif c==delimiter_n
                push!(row, "")
                # state remains the same
            else
            	n_field += 1
                field[n_field] = c
                state = FieldState.in_un_quoted_field
            end
        elseif state==FieldState.in_quoted_field
            if c==quotechar_n
                state = FieldState.in_double_quoted_field
            else
                if n_field==length(field)
                    resize!(field, length(field)*2)
                end
                n_field += 1
                field[n_field] = c # state remains the same
            end
        elseif state==FieldState.in_un_quoted_field
            if c==r_newline_n || c==n_newline_n
            	resize!(field, n_field)
                push!(row, string_type(field))
                produce(row)
                field = Array(Uint8, 1024)
                n_field = 0
                row = Array(string_type, 0)
                state = FieldState.not_in_field
            elseif c==delimiter_n
                resize!(field, n_field)
                push!(row, string_type(field))
                field = Array(Uint8, 1024)
                n_field = 0
                state = FieldState.not_in_field
            else
                if n_field==length(field)
                    resize!(field, length(field)*2)
                end
                n_field += 1
                field[n_field] = c # state remains the same
            end
        elseif state==FieldState.in_double_quoted_field
            if c==quotechar_n
                if n_field==length(field)
                    resize!(field, length(field)*2)
                end
                n_field += 1
                field[n_field] = c
                state = FieldState.in_quoted_field
            elseif c==r_newline_n || c==n_newline_n
                resize!(field, n_field)
                push!(row, string_type(field))
                field = Array(Uint8, 1024)
                n_field = 0
                produce(row)
                row = Array(string_type, 0)
                state = FieldState.not_in_field
            elseif c==delimiter_n
                resize!(field, n_field)
                push!(row, string_type(field))
                field = Array(Uint8, 1024)
                n_field = 0
                state = FieldState.not_in_field
            else
                throw(@sprintf("Invalid Byte In Double Quoted Field: %s", string(c)))
            end
        else
            throw("Invalid State")
        end
    end
    if n_field>0
        resize!(field, n_field)
        push!(row, string_type(field))
    end
    if length(row)>0
    	produce(row)
    end
    if close_io
    	close(io)
    end
end

function clean_col_names(field_names::Array{ASCIIString,1})
    cleaned_field_names = Array(ASCIIString, 0)
    for field_name=field_names
        cleaned_field_name = field_name
        if "type"==field_name
            cleaned_field_name = "type_0"
        end
        push!(cleaned_field_names, cleaned_field_name)
    end
    cleaned_field_names
end

abstract XsvRow
function create_xsv_row_type(field_names::Array{ASCIIString,1}; strtype="ascii"::ASCIIString, type_name=""::ASCIIString)
    field_names = clean_col_names(field_names)
    string_type = get_string_type(strtype::String)
    if type_name==""
        num = 1
        while isdefined(symbol("Xsv" * string(num) * "Row"))
            num += 1
        end
        type_name = "Xsv" * string(num) * "Row"
    end
    type_string = "type $type_name <: XsvRow\n";

    for field_name = field_names
        type_string = type_string*"    $field_name::$string_type\n"
    end
    eval(parse(type_string*"end"))
    eval(parse(type_name))
end

iterxsv(io::IO; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String, close_io=false::Bool) = @task _iterxsv(io, delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=close_io)
iterxsv(data::String; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String) = iterxsv(IOBuffer(data), delimiter=delimiter, quotechar=quotechar, strtype=strtype)
iterxsv(data::String; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String) = iterxsv(IOBuffer(data), delimiter=delimiter, quotechar=quotechar, strtype=strtype)
function _iterxsvh(io::IO; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String, close_io=false::Bool)
     xsv_stream = iterxsv(io, delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=close_io)
     header = consume(xsv_stream)
     for row in xsv_stream
        produce([header[i]=>row[i] for i=1:length(header)])
     end
end
iterxsvh(data::String; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String) = @task _iterxsvh(IOBuffer(data), delimiter=delimiter, quotechar=quotechar, strtype=strtype)
function _iterxsvht(io::IO; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String, close_io=false::Bool, row_type_name=""::ASCIIString)
    xsv_stream = iterxsv(io, delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=close_io)
    header = consume(xsv_stream)
    xsv_row = create_xsv_row_type(header, strtype=strtype, type_name=row_type_name)
    for row in xsv_stream
       produce(xsv_row(row...))
    end
end
iterxsvht(data::String; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String, row_type_name=""::ASCIIString) = @task _iterxsvht(IOBuffer(data), delimiter=delimiter, quotechar=quotechar, strtype=strtype, row_type_name=row_type_name)

function readxsv(io::IO; delimiter=','::Char, quotechar='"'::Char, strtype="ascii"::String, close_io=false::Bool)
	string_type = get_string_type(strtype)
	rows = Array(Vector{string_type}, 0)
	for row in iterxsv(io, delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=close_io)
		push!(rows, row)
	end
	rows
end
readxsv(data::String; delimiter=',', quotechar='"', strtype="ascii"::String) = readxsv(IOBuffer(data), delimiter=delimiter, quotechar=quotechar, strtype=strtype)

fiterxsv(xsv_file::String; delimiter=',', quotechar='"', strtype="ascii"::String) = @task _iterxsv(open(xsv_file), delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=true)
fiterxsvh(xsv_file::String; delimiter=',', quotechar='"', strtype="ascii"::String) = @task _iterxsvh(open(xsv_file), delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=true)
fiterxsvht(xsv_file::String; delimiter=',', quotechar='"', strtype="ascii"::String) = @task _iterxsvht(open(xsv_file), delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=true)
freadxsv(xsv_file::String; delimiter=',', quotechar='"', strtype="ascii"::String) = readxsv(open(xsv_file), delimiter=delimiter, quotechar=quotechar, strtype=strtype, close_io=true)


text:
# ProgressBarLogger

"""
    ProgressBarLogger()

Prints the progress of the current epoch.
"""
mutable struct ProgressBarLogger <: SafeCallback
    p::Union{Nothing,Progress}
end
ProgressBarLogger() = ProgressBarLogger(nothing)

function on(::EpochBegin,
        phase::AbstractFittingPhase,
        cb::ProgressBarLogger,
        learner)
    e = learner.cbstate[:history].epochs + 1
    cb.p = Progress(numsteps(learner, phase), "Epoch $(e) $(phase): ")
end

on(::BatchEnd, ::AbstractFittingPhase, cb::ProgressBarLogger, learner) = next!(cb.p)

runafter(::ProgressBarLogger) = (Recorder,)
stateaccess(::ProgressBarLogger) = (data = Read(), cbstate = (history = Read()),)

# MetricsLogger

"""
    MetricsLogger()

Prints the metrics after every epoch.
"""
struct MetricsLogger <: SafeCallback end

function on(::EpochEnd,
        phase::AbstractFittingPhase,
        cb::MetricsLogger,
        learner)

    for metric in getmetrics(learner.callbacks)
        println(string(metric), ": ", epochvalue(metric))
    end
end

stateaccess(::MetricsLogger) = (callbacks = Read(),)
runafter(::MetricsLogger) = (AbstractMetric,)

# StopOnNaNLoss

"""
    StopOnNaNLoss()

Stops the training when a NaN loss is encountered.
"""
struct StopOnNaNLoss <: SafeCallback end

function on(::BackwardEnd, ::AbstractTrainingPhase, ::StopOnNaNLoss, learner)
    !isnan(learner.batch.loss) || throw(CancelFittingException("Encountered NaN loss"))
end

stateaccess(::StopOnNaNLoss) = (batch = (loss = Read()),)

# Early stopping
mutable struct EarlyStopping <: SafeCallback
    patience::Int
    waited::Int
    lowest::Float64
end
EarlyStopping(patience) = EarlyStopping(patience, 0, Inf64)

function on(::EpochEnd, ::ValidationPhase, cb::EarlyStopping, learner)
    valloss = epochvalue(getloss(learner.callbacks))
    if (valloss > cb.lowest)
        if !(cb.waited < cb.patience)
            throw(CancelFittingException("Validation loss did not improve for $(cb.patience) epochs"))
        else
            cb.waited += 1
        end
    else
        cb.waited = 0
        cb.lowest = valloss
    end
end

stateaccess(::EarlyStopping) = (callbacks = Read(),)

struct ToGPU <: SafeCallback end

function on(::EpochBegin, ::AbstractFittingPhase, ::ToGPU, learner)
    model!(learner, gpu(learner.model))
end

stateaccess(::ToGPU) = (model = Write(), params = Write(), batch = (xs = Write(), ys = Write()),)

function on(::BatchBegin, ::AbstractFittingPhase, cb::ToGPU, learner)
    learner.batch.xs = gpu(learner.batch.xs)
    learner.batch.ys = gpu(learner.batch.ys)
end

garbagecollect() = (GC.gc(); ccall(:malloc_trim, Cvoid, (Cint,), 0))

"""
    GarbageCollect(nsteps)

Every `nsteps` steps, forces garbage collection.
Use this if you get memory leaks from, for example, parallel data loading.
"""
function GarbageCollect(nsteps::Int = 100)
    return CustomCallback{BatchEnd, AbstractFittingPhase}(nsteps) do learner
        garbagecollect()
    end
end

ext: jl


text:
include("../include.jl")
include("create_report.jl")
overlapping_results = get_CUTEst_results()

####
#### Lets do an overall summary of the results
####

# i.e. create Table 3 through 5.

# compute median number of iterations

##### Comparison on time to find a KKT point
#####
##### Lets take a look at the problems for which
##### - both solvers return a KKT point with the same objective
#####
#####
##### - one of the solvers finds a KKT point with the same objective value
##### - neither finds an unbounded solution
#####

f_TOL = 1e-1

function create_opt_res(res::Dict{String, Dict{String,problem_summary}}, method_list)
    opt_res = Dict{String, Dict{String,problem_summary}}()
    for method in method_list
        opt_res[method] = Dict{String,problem_summary}()
    end

    problem_list = keys(first(res)[2])
    for problem_name in problem_list
      same_fval = true
      first = true
      fval_others = Inf
      some_method_failed = false

      for method in method_list
         fval = res[method][problem_name].fval
         if first || abs(fval_others - fval) / (1.0 + max(abs(fval),abs(fval_others))) < f_TOL
           fval_others = min(fval_others, fval)
           first = false
         else
           same_fval = false
         end
         if res[method][problem_name].status != :Optimal
           some_method_failed = true
         end
      end

      if same_fval || some_method_failed
          for method in method_list
             opt_res[method][problem_name] = res[method][problem_name]
          end
      end
    end
    return opt_res
end

method_list = keys(overlapping_opt_results)
opt_res = create_opt_res(overlapping_opt_results, method_list)

its, best, ratios, times = compute_its_etc(opt_res,MAX_IT=3000);

using PyPlot

plot_iteration_ratios(its, best, ratios)
savefig("output/opt_iter_ratios.pdf")

plot_iterations(its, best, ratios, 3000)
savefig("output/opt_iter_curve.pdf")

##
## FUNCTION VALUES
##

function best_fvals(method_name1, method_name2, overlapping_opt_results)
    better_fval = Dict()
    succeed = Dict()

    feas_tol = 1e-4
    fval_tol = 1e-1

    better_fval[method_name1] = 0
    better_fval[method_name2] = 0

    succeed[method_name1] = 0
    succeed[method_name2] = 0

    problem_list = keys(first(overlapping_opt_results)[2])

    for problem_name in problem_list
      fval1 = overlapping_opt_results[method_name1][problem_name].fval
      con_vio1 = overlapping_opt_results[method_name1][problem_name].con_vio

      fval2 = overlapping_opt_results[method_name2][problem_name].fval
      con_vio2 = overlapping_opt_results[method_name2][problem_name].con_vio

      if con_vio1 < feas_tol && con_vio2 < feas_tol
        fbest = min(fval1,fval2)
        ftol_scaled = fval_tol * (1.0 + abs(fbest))

        if fval1 < fval2 - ftol_scaled
          better_fval[method_name1] += 1
        elseif fval2 < fval1 - ftol_scaled
          better_fval[method_name2] += 1
        end
      end

      if con_vio1 < feas_tol
        succeed[method_name1] += 1
      end

      if con_vio2 < feas_tol
        succeed[method_name2] += 1
      end
    end

    return better_fval, succeed
end

method_name1 = "ipopt"
method_name2 = "one phase"
best_fvals(method_name1, method_name2, overlapping_opt_results)

##
## PLOT TIME CURVES
##
best = best_its(times)

min_y = 1.0

ratios = Dict()
y_vals = collect(1:length(best)) / length(best)
for (method_name, val) in its
  ratios[method_name] = times[method_name] ./ best
  semilogx(sort(ratios[method_name]), y_vals, label=method_name, basex=2)
  min_y = min(min_y,sum(ratios[method_name] .== 1.0) / length(best))
end
ax = gca()
ax[:set_xlim]([1.0,2^8.0])
ax[:set_ylim]([min_y,1.0])

#ax[:xaxis][:ticker] = 0.5
#title("Comparsion on 45 CUTEst problems")
xlabel("times ratio to best solver")
ylabel("proportion of problems")

legend()

function create_opt_res(res::Dict{String, Dict{String,problem_summary}}, method_list; f_TOL = 1e-1)
    opt_res = Dict{String, Dict{String,problem_summary}}()
    for method in method_list
        opt_res[method] = Dict{String,problem_summary}()
    end

    problem_list = keys(first(res)[2])
    for problem_name in problem_list
      same_fval = true
      first = true
      fval_others = Inf

      for method in method_list
         fval = res[method][problem_name].fval
         if res[method][problem_name].status == :Optimal
           if first || abs(fval_others - fval) / (1.0 + max(abs(fval),abs(fval_others))) < f_TOL
             fval_others = min(fval_others, fval)
             first = false
           else
             same_fval = false
           end
         else
           same_fval = false
         end
      end

      if same_fval #|| some_method_failed
          for method in method_list
             opt_res[method][problem_name] = res[method][problem_name]
          end
      end
    end
    return opt_res
end

method_list = keys(overlapping_opt_results)
#opt_res = create_opt_res(overlapping_opt_results, method_list, f_TOL = 1e-2)


text:
function testTrain(train, topol; toroidal = false,
                    normType = :zscore, kernel = gaussianKernel)

    xdim = 8
    ydim = 10

    som = initSOM(train, xdim, ydim, norm = normType,
                  topol = topol, toroidal = toroidal)
    som = trainSOM(som, train, 10000, kernelFun = kernel)

    ntrain = nrow(train)
    npopul = sum(som.population)

    return ntrain == npopul
end

function testVisual(train, topol)

    xdim = 8
    ydim = 10

    tr  = train[1:50,:]
    val = train[51:100,:]

    som = initSOM(tr, xdim, ydim, norm = :zscore,
    topol = topol, toroidal = false)
    som = trainSOM(som, tr, 10000)

    vis = mapToSOM(som, tr)

    ntr    = nrow(tr)
    npopul = nrow(vis)

    return ntr == npopul
end

function testFreqs(train, wClasses, classes)

    xdim = 8
    ydim = 10

    som = initSOM(train, xdim, ydim)
    som = trainSOM(som, train, 10000)

    f = classFrequencies(som, wClasses, classes)

    res = sort( [String(x) for x in names(f)[5:end]])
    ori = sort(unique(iris[:Species]))

    return res == ori
end

function testDensityPlot(train, topol)

    xdim = 8
    ydim = 10
    fName = "density.png"

    som = initSOM(train, xdim, ydim, topol = topol)
    som = trainSOM(som, train, 10000)

    if MPL_INSTALLED
        plotDensity(som, fileName = fName)
        ok = isfile(fName)
    else
        result = plotDensity(som, fileName = fName)
        ok = typeof(result) == Symbol && result == :ERR_MPL
    end

    return ok
end

function testOtherDensityPlot(tr, pred)

    xdim = 8
    ydim = 10
    fName = "density.png"

    som = initSOM(tr, xdim, ydim)
    som = trainSOM(som, tr, 10000)

    vis = mapToSOM(som, pred)

    if MPL_INSTALLED
        plotDensity(som, predict = vis, fileName = fName)
        ok = isfile(fName)
    else
        result = plotDensity(som, predict = vis, fileName = fName)

        ok = typeof(result) == Symbol && result == :ERR_MPL
    end

    return ok
end

function testClassesPlot(train, wClasses, topol)

    xdim = 8
    ydim = 10
    fName = "classes.png"

    som = initSOM(train, xdim, ydim, topol = topol)
    som = trainSOM(som, train, 10000)

    freqs = classFrequencies(som, wClasses, :Species)
    if MPL_INSTALLED
        plotClasses(som, freqs, fileName = fName)
        ok = isfile(fName)
    else
        result = plotClasses(som, freqs, fileName = fName)

        ok = typeof(result) == Symbol && result == :ERR_MPL
    end

    return ok
end

function testClassesToFile(train, wClasses)

    xdim = 8
    ydim = 10
    fName = "classesfile"

    som = initSOM(train, xdim, ydim)
    som = trainSOM(som, train, 10000)

    freqs = classFrequencies(som, wClasses, :Species)
    if MPL_INSTALLED
        plotClasses(som, freqs, device = :svg, fileName = fName)
        ok = isfile(fName * ".svg")
    else
        result = plotClasses(som, freqs, device = :svg, fileName = fName)

        ok = typeof(result) == Symbol && result == :ERR_MPL
    end

    return ok
end


text:
using NCDatasets
using Test
using DataStructures
using Dates

time = DateTime(2000,1,1):Dates.Day(1):DateTime(2000,1,3)
time_bounds = Matrix{DateTime}(undef,2,length(time))
time_bounds[1,:] = time .- Dates.Hour(12)
time_bounds[2,:] = time .+ Dates.Hour(12)

fname = tempname()
ds = NCDataset(fname,"c")
ds.dim["nv"] = 2;
ds.dim["time"] = length(time);

nctime = defVar(ds, "time", Float64, ("time",),attrib=OrderedDict(
    "units" => "days since 2000-01-01",
    "scale_factor" => 10.,
    "bounds" => "time_bounds"));

nctime_bounds = defVar(ds, "time_bounds", Float64, ("nv","time"),attrib=OrderedDict())

ncvar = nctime
nctime_bounds = NCDatasets.bounds(ncvar)
nctime_bounds[:,:] = time_bounds
@test nctime_bounds.var[:,:] ≈ [-0.5  0.5  1.5;
                                0.5  1.5  2.5]
@test nctime_bounds[:,:] == time_bounds

close(ds)


text:
x = 1; y = 2

if (x < y)
  println("x($x) is less than y($y)")
elseif (x > y)
  println("x($x) is greater than y($y)")
else
  println("x($x) is equal to y($y)")
end

# Leaky Variable Scope

if (x < y)
  relation = "x($x) is less than y($y)"
elseif (x > y)
  relation = "x($x) is greater than y($y)"
else
  relation = "x($x) is equal to y($y)"
end

println(relation)

# Ternary Operator
# - can be nested

z = x == y ? "x($x) is equal to y($y)" : "x($x) is not equal to y($y)"

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
using ApproxFun, SingularIntegralEquations, Base.Test

println("Hilbert test")

x=Fun(identity)
f=(exp(x)/(sqrt(1-x)*sqrt(x+1)))
@test_approx_eq (Hilbert(f|>space,0)*f)(.1) (-0.8545003781055088)
@test_approx_eq (Hilbert(0)*f)(.1) (-0.8545003781055088)
@test_approx_eq (Hilbert()*f)(.1) 1.1404096104609646386

x=Fun(identity,[-1,2])
f=(exp(x)/(sqrt(2-x)*sqrt(x+1)))
@test_approx_eq (Hilbert(f|>space,0)*f)(.1) 0.49127801561694168644
@test_approx_eq (Hilbert(0)*f)(.1) 0.49127801561694168644
@test_approx_eq (Hilbert()*f)(.1) 1.6649936695644078289

x=Fun(identity)
f=(exp(x)*(sqrt(1-x)*sqrt(x+1)))
@test_approx_eq (Hilbert()*f)(.1) 0.43723982258866913063

x=Fun(identity,[-1,2])
f=(exp(x)*(sqrt(2-x)*sqrt(x+1)))
@test_approx_eq (Hilbert()*f)(.1) 2.1380903070701673244

x=Fun(identity)
d=domain(x)
w=1/sqrt(1-x^2)
H=Hilbert(d)
B=ldirichlet(d)

for a in [sqrt(sqrt(5)-2)/2,1.,10.]
    L=H[w]+1/a/sqrt(1+a^2)*x
    u=[B,L]\[1.]
    usol = (1+a^2)/(x^2+a^2)
    @test norm(u-usol) <= eps(1000/a)
end

x = Fun(identity)
w = 1/sqrt(1-x^2)
H = Hilbert(space(w))
@test_approx_eq  (H[w]*exp(x))(.1) hilbert(w*exp(x))(.1)

x = Fun(identity)
w = sqrt(1-x^2)
H = Hilbert(space(w))
@test_approx_eq (H[w]*exp(x))(.1) hilbert(w*exp(x))(.1)

println("Stieltjes test")

ds1 = JacobiWeight(-.5,-.5,ApproxFun.ChebyshevDirichlet{1,1}())
ds2 = JacobiWeight(-.5,-.5,Chebyshev())
rs = Chebyshev([2.,4.+3im])
f1 = Fun(x->exp(x)/sqrt(1-x^2),ds1)
f2 = Fun(x->exp(x)/sqrt(1-x^2),ds2)
S = Stieltjes(ds1,rs)

z = 3.+1.5im
@test_approx_eq (S*f1)(z) stieltjes(f2,z) #val,err = quadgk(x->f1(x)./(z-x),-1.,1.)
# Operator 1.1589646343327578 - 0.7273679005911196im
# Function 1.1589646343327455 - 0.7273679005911283im

ds1 = JacobiWeight(.5,.5,Ultraspherical(1))
ds2 = JacobiWeight(.5,.5,Chebyshev())
rs = Chebyshev([2.,4.])
f1 = Fun(x->exp(x)*sqrt(1-x^2),ds1)
f2 = Fun(x->exp(x)*sqrt(1-x^2),ds2)
S = Stieltjes(ds1,rs)

z = 3.
@test_approx_eq (S*f1)(z) stieltjes(f2,z) #val,err = quadgk(x->f1(x)./(z-x),-1.,1.;reltol=eps())
# Operator 0.6616422557285478 + 0.0im
# Function 0.661642255728541 - 0.0im

println("Stieltjes integral test")

ds1 = JacobiWeight(-.5,-.5,ApproxFun.ChebyshevDirichlet{1,1}())
ds2 = JacobiWeight(-.5,-.5,Chebyshev())
rs = Chebyshev([2.,4.])
f1 = Fun(x->exp(x)/sqrt(1-x^2),ds1)
f2 = Fun(x->exp(x)/sqrt(1-x^2),ds2)
S = Stieltjes(ds1,rs,0)

z = 3.
@test_approx_eq (S*f1)(z) SingularIntegralEquations.stieltjesintegral(f2,z)
# Operator 3.6322473044237698 + 0.0im
# Function 3.6322473044237515

ds1 = JacobiWeight(.5,.5,Ultraspherical(1))
ds2 = JacobiWeight(.5,.5,Chebyshev())
rs = Chebyshev([2.0,4.0])
f1 = Fun(x->exp(x)*sqrt(1-x^2),ds1)
f2 = Fun(x->exp(x)*sqrt(1-x^2),ds2)
S = Stieltjes(ds1,rs,0)

z = 3.0
@test_approx_eq (S*f1)(z) SingularIntegralEquations.stieltjesintegral(f2,z)
# Operator 1.7772163062194861 + 0.0im
# Function 1.7772163062194637

println("Cauchy test")

f2=Fun(sech,PeriodicLine())
@test_approx_eq cauchy(f2,1.+im) (0.23294739894134472 + 0.10998776661109881im )
@test_approx_eq cauchy(f2,1.-im) (-0.23294739894134472 + 0.10998776661109881im )

f=Fun(sech,Line())
@test_approx_eq cauchy(f,1.+im) cauchy(f2,1.+im)

f=Fun(z->exp(exp(0.1im)*z+1/z),Laurent(Circle()))

@test_approx_eq hilbert(f,exp(0.2im)) hilbert(Fun(f,Fourier),exp(0.2im))
@test_approx_eq hilbert(f,exp(0.2im)) -hilbert(reverseorientation(f),exp(0.2im))
@test_approx_eq hilbert(f,exp(0.2im)) -hilbert(reverseorientation(Fun(f,Fourier)),exp(0.2im))

@test_approx_eq cauchy(f,0.5exp(0.2im)) cauchy(Fun(f,Fourier),0.5exp(0.2im))
@test_approx_eq cauchy(f,0.5exp(0.2im)) -cauchy(reverseorientation(f),0.5exp(0.2im))
@test_approx_eq cauchy(f,0.5exp(0.2im)) -cauchy(reverseorientation(Fun(f,Fourier)),0.5exp(0.2im))
@test_approx_eq cauchy(f,0.5exp(0.2im)) (OffHilbert(space(f),Laurent(Circle(0.5)))*f)(0.5exp(0.2im))/(2im)

f=Fun(z->exp(exp(0.1im)*z+1/(z-1.)),Laurent(Circle(1.,0.5)))
@test_approx_eq cauchy(f,0.5exp(0.2im)) cauchy(Fun(f,Fourier),0.5exp(0.2im))
@test_approx_eq cauchy(f,0.5exp(0.2im)) -cauchy(reverseorientation(f),0.5exp(0.2im))
@test_approx_eq cauchy(f,0.5exp(0.2im)) -cauchy(reverseorientation(Fun(f,Fourier)),0.5exp(0.2im))

Γ=Circle()∪Circle(0.5)
f=depiece([Fun(z->z^(-1),Γ[1]),Fun(z->z,Γ[2])])
A=I-(f-Fun(one,space(f)))*Cauchy(-1)
u=A\(f-Fun(one,space(f)))

@test_approx_eq 1+cauchy(u,.1) 1
@test_approx_eq 1+cauchy(u,.8) 1/0.8
@test_approx_eq 1+cauchy(u,2.) 1

c1=0.5+0.1;r1=3.;
c2=-0.1+.2im;r2=0.3;
d1=Circle(c1,r1)
d2=Circle(c2,r2)
z=Fun(identity,d2);
C=Cauchy(Space(d1),Space(d2))
@test norm((C*Fun(exp,d1)-Fun(exp,d2)).coefficients)<100eps()

C2=Cauchy(Space(d2),Space(d1))
@test norm((C2*Fun(z->exp(1/z)-1,d2)+Fun(z->exp(1/z)-1,d1)).coefficients)<100000eps()

c1=0.1+.1im;r1=.4;
c2=-2.+.2im;r2=0.3;
d1=Circle(c1,r1)
d2=Circle(c2,r2)
@test norm((Cauchy(d1,d2)*Fun(z->exp(1/z)-1,d1)+Fun(z->exp(1/z)-1,d2)).coefficients)<2000eps()

# complex contour

#Legendre uses FastGuassQuadrature
f=Fun(exp,Legendre())
#@test_approx_eq cauchy(f,.1+0.000000000001im) cauchy(f,.1,+)
#@test_approx_eq cauchy(f,.1-0.000000000001im) cauchy(f,.1,-)
#@test_approx_eq (cauchy(f,.1,+)-cauchy(f,.1,-)) exp(.1)

ω=2.
d=Interval(0.5im,30.0im/ω)
x=Fun(identity,Legendre(d))
@test_approx_eq cauchy(exp(im*ω*x),1.+im) (-0.025430235512791915911 + 0.0016246822285867573678im)

println("Arc test")

a=Arc(0.,1.,0.,π/2)
ζ=Fun(identity,a)
f=Fun(exp,a)*sqrt(abs((ζ-1)*(ζ-im)))
z=.1+.2im
#@test_approx_eq cauchy(f,z) sum(f/(ζ-z))/(2π*im)
z=exp(.1im)
@test_approx_eq hilbert(f,z) im*(cauchy(f,z,+)+cauchy(f,z,-))

println("Functional test")
z=.1+.2im
x=Fun(identity)
f=exp(x)*sqrt(1-x^2)
@test_approx_eq Stieltjes(space(f),z)*f stieltjes(f,z)

a=Arc(0.,1.,0.,π/2)
ζ=Fun(identity,a)
f=Fun(exp,a)*sqrt(abs((ζ-1)*(ζ-im)))
H=Hilbert()
z=exp(.1im)
@test_approx_eq (H*f)(z) hilbert(f,z)

println("Logkernel test")

a=1.0;b=2.0
d=Interval(a,b)
z=Fun(d)
f=real(exp(z)/(sqrt(z-a)*sqrt(b-z)))
S=space(f)
x=4.0+2im;
@test_approx_eq linesum(f*log(abs(x-z))) logkernel(f,x)*π

a=1.0;b=2.0+im
d=Interval(a,b)
z=Fun(d)
f=real(exp(z)/(sqrt(z-a)*sqrt(b-z)))
S=space(f)
x=4.0+2im;
@test_approx_eq linesum(f*log(abs(x-z))) 13.740676344264614
@test_approx_eq linesum(f*log(abs(x-z))) logkernel(f,x)*π

linesum(f)
sum(f)
a=1.0;b=2.0+im
d=Interval(a,b)
z=Fun(d)
f=real(exp(z)*(sqrt(z-a)*sqrt(b-z)))
S=space(f)
x=4.0+2im;
@test_approx_eq linesum(f*log(abs(x-z))) logkernel(f,x)*π

a=1.0;b=2.0
d=Interval(a,b)
z=Fun(d)
f=real(exp(z)/(sqrt(z-a)*sqrt(b-z)))
x=1.5
@test_approx_eq (SingularIntegral(space(f),0)*f)(x) logkernel(f,x)

f=real(exp(z)*(sqrt(z-a)*sqrt(b-z)))
x=1.5
@test_approx_eq (SingularIntegral(space(f),0)*f)(x) logkernel(f,x)

a=1.0;b=2.0+im
d=Interval(a,b)
z=Fun(d)
f=real(exp(z)/(sqrt(z-a)*sqrt(b-z)))
x=1.5+0.5im
H=SingularIntegral(space(f),0)
@test_approx_eq (H*f)(x) logkernel(f,x)

f=real(exp(z)*(sqrt(z-a)*sqrt(b-z)))
@test_approx_eq (SingularIntegral(space(f),0)*f)(x) logkernel(f,x)

a=1.0;b=2.0+im
d=Interval(a,b)
z=Fun(d)
f=real(exp(z)/(sqrt(z-a)*sqrt(b-z)))
S=JacobiWeight(-0.5,-0.5,ChebyshevDirichlet{1,1}(d))
H=OffSingularIntegral(S,Chebyshev([3,4]),0)
@test_approx_eq (H*f)(3.5) logkernel(f,3.5)

H=OffSingularIntegral(S,Chebyshev([3,4.0+im]),0)
@test_approx_eq (H*f)(3.5+0.5im) logkernel(f,3.5+0.5im)

## Circle

d=Circle(0.2,3.0)
S=Fourier(d)
ζ=Fun(d)
f=real(ζ+1/(ζ-0.1))
#z=0.1+0.1im;@test_approx_eq linesum(log(abs(ζ-z))*f) logkernel(f,z)*π
#z=5.0+0.1im;@test_approx_eq linesum(log(abs(ζ-z))*f) logkernel(f,z)*π

d=Circle(0.2,3.0)
S=Fourier(d)
H=Hilbert(S,0)
ζ=Fun(d)
f=real(ζ+1/(ζ-0.1))
z=0.2+3im;@test_approx_eq (H*f)(z) logkernel(f,z)

if isdir(Pkg.dir("FastGaussQuadrature"))
    println("Stieltjes moment test")
    include("stieltjesmomentTest.jl")
end

println("Curve Test")

include("CurveTest.jl")

println("FundamentalSolutions test")

include("FundamentalSolutionsTest.jl")

println("Convolution ProductFun test")

include("convolutionProductFunTest.jl")

println("LowRankMatrix test")

include("LowRankMatrixTest.jl")

println("HierarchicalVector test")

include("HierarchicalVectorTest.jl")

println("Hierarchical solve test")

include("hierarchicalsolveTest.jl")

# Testing Travis CI abilities. Currently, Examples tests are slow.
#println("\nExamples test\n")
#
#include("ExamplesTest.jl")

println("WienerHopfTest")
include("WienerHopfTest.jl")


text:
nb = 4

#initialisation variables
x_a = zeros(Float64 , nb) ; v_a = zeros(Float64 , nb) ; p_a = zeros(Float64 , nb)
f_a = zeros(Float64, nb)
x_b = zeros(Float64 , nb) ; v_b = zeros(Float64 , nb) ; p_b = zeros(Float64 , nb)
f_b = zeros(Float64, nb)
g = zeros(Float64 , nb)

#intiialisation valeurs
rng = [0.679923 0.269864 0.88061 0.722799 0.566828 0.44927 0.935181 0.1669 0.29829 0.670977 0.369223 0.927932]
x_a[1] = 0.5 ; v_a[1] = 0
x_b[1] = 1.25 ; v_b[1] = 0
w = 1 ; c1 = 2 ; c2 = 2

function swarm(x_a,v_a,p_a,x_b,v_b,p_b,g,rng,w,c1,c2)
    cpt = 1
    for t in 1:nb
        push!(f_b, round(f(x_b[t]); digits = 4))
        if t != 1
            v_a[t] = velocite(w, v_a[t-1], rng[cpt], c1, p_a[t-1], x_a[t-1], rng[cpt + 1], c2, g[t-1])
            cpt += 2
            v_b[t] = velocite(w, v_b[t-1], rng[cpt], c1, p_b[t-1], x_b[t-1], rng[cpt + 1], c2, g[t-1])
            cpt += 2
            x_a[t] = position(x_a[t-1], v_a[t])
            x_b[t] = position(x_b[t-1], v_b[t])
        end
        f_a[t] = f(round(x_a[t];digits=4))
        f_b[t] = f(round(x_b[t];digits=4))
        p_a[t] = perso_max(x_a)
        p_b[t] = perso_max(x_b)
        g[t] = pos_max(x_a[t],x_b[t])
    end
    println("x_a: ", x_a) ; println("v_a: ", v_a) ; println("p_a: ", p_a)
    println("f_a: ", f_a)
    println("g: ", g)
    println("x_b: ", x_b); println("v_b: ", v_b); println("p_b: ", p_b)
    println("f_b: ", f_b)
end

function f(x)
    return(x*sin(pi*10*x)+1)
end

function velocite(w, v, rng1, c1, p, x, rng2, c2, g)
    return((w*v) + (rng1*c1*(p - x)) + (rng2*c2*(g - x)))
end

function position(x,v)
    return(x+v)
end

function pos_max(x1,x2)
    if f(x1)-f(x2) > 0
        return(x1)
    else
        return(x2)
    end
end

function perso_max(x)
    max = f(x[1])
    pos = 1
    for i in 2:length(x)
        if f(x[i]) > max
            max = f(x[i])
            pos = i
        end
    end
    return(x[pos])
end

swarm(x_a, v_a, p_a, x_b, v_b, p_b, g, rng, w, c1, c2)


text:
########################################################################################
# timed pmapsum and pmapreduce
########################################################################################
function pmapsum(comm::MPI.Comm, f, iters, args...; kwargs...)
	print_timings = get(kwargs, :print_timings, false)
	localtimer = TimerOutput()
	ℓ_ωind_iter_on_proc = productsplit(iters, comm)
	res_local = f(localtimer, ℓ_ωind_iter_on_proc, args...; kwargs...)
	if print_timings && res_local !== nothing
		println(localtimer)
	end
	MPI.Reduce(res_local, MPI.SUM, 0, comm)
end

function pmapsum(::Any, f, iters, args...; kwargs...)
	print_timings = get(kwargs, :print_timings, false)
	fwrapper(ps) = begin
		localtimer = TimerOutput()
		res = f(localtimer, ps, args...)
		if print_timings && ParallelUtilities.workerrank(ps) == 1
			println(localtimer)
		end
		res
	end

	pmapreduce_productsplit(fwrapper, ParallelUtilities.elementwisesum!, iters...)
end

function pmapreduce(::Any, f, fred, iters, args...;	kwargs...)
	print_timings = get(kwargs, :print_timings, false)
	fwrapper(ps) = begin
		localtimer = TimerOutput()
		res = f(localtimer, ps, args...)
		if print_timings && ParallelUtilities.workerrank(ps) == 1
			println(localtimer)
		end
		res
	end

	pmapreduce_productsplit(fwrapper, fred, iters...)
end

########################################################################################
# utility functions
########################################################################################

signaltomaster!(r::RemoteChannel{Channel{Bool}}) = put!(r, true)
signaltomaster!(r::RemoteChannel{Channel{TimerOutput}}, timer::TimerOutput) = put!(r, timer)
# Ignore anything other than RemoteChannels
signaltomaster!(val...) = nothing

_broadcast(A, root, comm) = A
_broadcast(A, root, comm::MPI.Comm) = MPI.bcast(A, root, comm)


text:
export buildPOP_1v1c, buildPOP_1v2c, buildPOP_1v2, buildPOP_EllJoszMolc
export buildPOP_WB2, buildPOP_WB5
export lasserre_ex1, lasserre_ex2, lasserre_ex3, lasserre_ex5
export get_WB5cliques, get_case9cliques

############################
### Geometric problems
############################

function buildPOP_1v1c()
    z = Variable("z", Complex)
    problem = Problem()
    add_variable!(problem, z)
    set_objective!(problem, -real(z))
    add_constraint!(problem, "ineq", abs2(z) << 4)
    return problem
end

function buildPOP_1v2c()
    z = Variable("z", Complex)
    problem = Problem()
    add_variable!(problem, z)
    set_objective!(problem, imag(z))
    add_constraint!(problem, "ineq_brn", abs2(z) << 1)
    θ = π/3
    add_constraint!(problem, "ineq_rot", real(z*exp(-im*θ)) >> 0)
    return problem
end

function buildPOP_1v2()
    x1 = Variable("x1", Real)
    x2 = Variable("x2", Real)
    problem = Problem()
    add_variable!(problem, x1); add_variable!(problem, x2)
    set_objective!(problem, -1.0*x1)
    add_constraint!(problem, "ineq", (x1^2+x2^2) << 1)
    θ1 = π/3
    add_constraint!(problem, "eq_rot1", (cos(θ1)*x1+sin(θ1)*x2) == 0)
    # θ2 = -π/3
    # add_constraint!(problem, "ineq_rot2", (cos(θ2)*x1+sin(θ2)*x2) >> 0)
    return problem
end

"""
    problem = buildPOP_2v3c

    Elliptic example problemp from Josz, Molzahn 2018 paper.
"""
function buildPOP_EllJoszMolc()
    z1 = Variable("z1", Complex)
    z2 = Variable("z2", Complex)
    problem = Problem()
    add_variable!(problem, z1); add_variable!(problem, z2);
    set_objective!(problem, 3-abs2(z1)-0.5*im*z1*conj(z2)^2+0.5im*z2^2*conj(z1))
    add_constraint!(problem, "eq1", (abs2(z1)-0.25*z1^2-0.25*conj(z1)^2) == 1)
    add_constraint!(problem, "eq2", (abs2(z1)+abs2(z2)) == 3)
    add_constraint!(problem, "eq3", (im*z2-im*conj(z2)) == 0)
    add_constraint!(problem, "ineq", (z2+conj(z2)) >> 0)
    return problem
end

############################
### OPF problems
############################

function buildPOP_WB2(; v2max = 0.976, rmeqs = false, setnetworkphase=false, addball=false)
    # OPFpbs = load_OPFproblems(MatpowerInput, joinpath("..", "data", "data_Matpower", "matpower", "WB2.m"))
    # problem_c = build_globalpb!(OPFpbs)
    problem_c, pt = import_from_dat(getinstancepath("Matpower", "QCQP", "WB2"))

    ## Converting to real ineq. only problem
    !rmeqs || change_eq_to_ineq!(problem_c)
    problem = pb_cplx2real(problem_c)

    if setnetworkphase
        ## Fixing volt phase of last bus to 0
        lastctr = problem.constraints["BaseCase_2_Volt_VOLTM_Re"]
        rm_constraint!(problem, "BaseCase_2_Volt_VOLTM_Re")

        ## Setting imag part to 0
        # pt = Point(SortedDict(Variable("BaseCase_2_VOLT_Im", Real)=>0.0), isdense=true)
        # infer_problem!(problem, pt)

        add_constraint!(problem, "BaseCase_2_Volt_VOLTM_Re", sqrt(lastctr.lb) << Variable("BaseCase_2_VOLT_Re", Real) << v2max)
        add_constraint!(problem, "BaseCase_2_Volt_VOLTM_Im", Variable("BaseCase_2_VOLT_Im", Real) == 0)
    elseif v2max != 0.976
        problem.constraints["BaseCase_2_Volt_VOLTM_Re"].ub = v2max^2
    end

    ## Adding ball constraint
    if addball
        p = Polynomial()
        for var in problem.variables
            p += Variable(var[1], var[2])^2
        end
        ub = problem.constraints["BaseCase_1_Volt_VOLTM_Re"].ub + v2max^2
        add_constraint!(problem, "Ball_ctr", p << ub)
    end

    return problem
end

function buildPOP_WB5(; q5min = 1.05, rmeqs = false)
    # OPFpbs = load_OPFproblems(MatpowerInput, joinpath("..", "data", "data_Matpower", "matpower", "WB5.m"))
    problem_c, pt = import_from_dat(getinstancepath("Matpower", "QCQP", "WB5"))

    # Sgen = OPFpbs["BaseCase"].ds.bus["BUS_5"]["Gen_1"].power_min
    # OPFpbs["BaseCase"].ds.bus["BUS_5"]["Gen_1"].power_min = real(Sgen) + im*q5min
    # problem_c = build_globalpb!(OPFpbs)

    ## Converting to real ineq. only problem
    !rmeqs || change_eq_to_ineq!(problem_c)
    return pb_cplx2real(problem_c)
end

############################
### Global Optim pbs from Lasserre2001
############################

"""
    problem, relax_ctx = lasserre_ex1()

    From Lasserre2001, global minimum : (3) -0.2428.
"""
function lasserre_ex1()
    x1 = Variable("x1", Real)
    x2 = Variable("x2", Real)
    problem = Problem()
    add_variable!(problem, x1); add_variable!(problem, x2)
    set_objective!(problem, (x1^2+1)^2 + (x2^2+1)^2 + (x1+x2+1)^2)

    relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                        d = 2)
    return problem, relax_ctx
end

"""
    problem, relax_ctx = lasserre_ex2()

    From Lasserre2001, global minimum : -11.4581.
"""
function lasserre_ex2()
    x1 = Variable("x1", Real)
    x2 = Variable("x2", Real)
    problem = Problem()
    add_variable!(problem, x1); add_variable!(problem, x2)
    set_objective!(problem, (x1^2+1)^2 + (x2^2+1)^2 -2*(x1+x2+1)^2)

    relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                        d = 2)
    return problem, relax_ctx
end

"""
    problem, relax_ctx = lasserre_ex3()

    From Lasserre2001, global minimum : -1/27, x1*² = x2*² = 1/3.
"""
function lasserre_ex3()
    x1 = Variable("x1", Real)
    x2 = Variable("x2", Real)
    problem = Problem()
    add_variable!(problem, x1); add_variable!(problem, x2)
    set_objective!(problem, x1^2 * x2^2 * (x1^2 + x2^2 - 1))

    relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                        d = 3)
    return problem, relax_ctx
end

"""
    problem, relax_ctx = lasserre_ex5()

    From Lasserre2001, global minimum : -2, for (1, 2).
    Relaxation : order 1 -> -3; order 2 -> -2.
"""
function lasserre_ex5(;d = 2)
    x1 = Variable("x1", Real)
    x2 = Variable("x2", Real)
    problem = Problem()
    add_variable!(problem, x1); add_variable!(problem, x2)
    set_objective!(problem, -(x1-1)^2 -(x1-x2)^2 -(x2-3)^2)
    add_constraint!(problem, "crt1", (1-(x1-1)^2) >> 0)
    add_constraint!(problem, "crt2", (1-(x1-x2)^2) >> 0)
    add_constraint!(problem, "crt3", (1-(x2-3)^2) >> 0)

    relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                        d = d)
    return problem, relax_ctx
end

function get_WB5cliques(relax_ctx, problem)
    if !relax_ctx.issparse
        return get_maxcliques(relax_ctx, problem)
    else
        maxcliques = Dict{String, Set{Variable}}()
        maxcliques["clique1"] = Set{Variable}([
            Variable("BaseCase_1_VOLT_Im", Real),
            Variable("BaseCase_1_VOLT_Re", Real),
            Variable("BaseCase_2_VOLT_Im", Real),
            Variable("BaseCase_2_VOLT_Re", Real),
            Variable("BaseCase_3_VOLT_Im", Real),
            Variable("BaseCase_3_VOLT_Re", Real)])
        maxcliques["clique2"] = Set{Variable}([
            Variable("BaseCase_2_VOLT_Im", Real),
            Variable("BaseCase_2_VOLT_Re", Real),
            Variable("BaseCase_3_VOLT_Im", Real),
            Variable("BaseCase_3_VOLT_Re", Real),
            Variable("BaseCase_4_VOLT_Im", Real),
            Variable("BaseCase_4_VOLT_Re", Real),
            Variable("BaseCase_5_VOLT_Im", Real),
            Variable("BaseCase_5_VOLT_Re", Real)])
        return maxcliques
    end
end

function get_case9cliques(relax_ctx, problem)
    if !relax_ctx.issparse
        return get_maxcliques(relax_ctx, problem)
    else
        maxcliques = Dict{String, Set{Variable}}()
        maxcliques["clique1"] = Set{Variable}([
            Variable("BaseCase_1_VOLT_Im", Real),
            Variable("BaseCase_1_VOLT_Re", Real),
            Variable("BaseCase_5_VOLT_Im", Real),
            Variable("BaseCase_5_VOLT_Re", Real),
            Variable("BaseCase_4_VOLT_Im", Real),
            Variable("BaseCase_4_VOLT_Re", Real),
            Variable("BaseCase_9_VOLT_Im", Real),
            Variable("BaseCase_9_VOLT_Re", Real),
            Variable("BaseCase_8_VOLT_Im", Real),
            Variable("BaseCase_8_VOLT_Re", Real)])
        maxcliques["clique2"] = Set{Variable}([
            Variable("BaseCase_2_VOLT_Im", Real),
            Variable("BaseCase_2_VOLT_Re", Real),
            Variable("BaseCase_9_VOLT_Im", Real),
            Variable("BaseCase_9_VOLT_Re", Real),
            Variable("BaseCase_8_VOLT_Im", Real),
            Variable("BaseCase_8_VOLT_Re", Real),
            Variable("BaseCase_7_VOLT_Im", Real),
            Variable("BaseCase_7_VOLT_Re", Real),
            Variable("BaseCase_6_VOLT_Im", Real),
            Variable("BaseCase_6_VOLT_Re", Real)])
        maxcliques["clique3"] = Set{Variable}([
            Variable("BaseCase_3_VOLT_Im", Real),
            Variable("BaseCase_3_VOLT_Re", Real),
            Variable("BaseCase_7_VOLT_Im", Real),
            Variable("BaseCase_7_VOLT_Re", Real),
            Variable("BaseCase_6_VOLT_Im", Real),
            Variable("BaseCase_6_VOLT_Re", Real),
            Variable("BaseCase_5_VOLT_Im", Real),
            Variable("BaseCase_5_VOLT_Re", Real),
            Variable("BaseCase_4_VOLT_Im", Real),
            Variable("BaseCase_4_VOLT_Re", Real)])
        return maxcliques
    end
end


text:
using Bridge, StaticArrays, Distributions
using Test, LinearAlgebra, Random

const percentile = 3.0

# Define a diffusion process
if ! @isdefined(OrnsteinUhlenbeck)
struct OrnsteinUhlenbeck  <: ContinuousTimeProcess{Float64}
    β::Float64 # drift parameter (also known as inverse relaxation time)
    σ::Float64 # diffusion parameter
    function OrnsteinUhlenbeck(β::Float64, σ::Float64)
        isnan(β) || β > 0.0 || error("Parameter λ must be positive.")
        isnan(σ) || σ > 0.0 || error("Parameter σ must be positive.")
        new(β, σ)
    end
end
end

# define drift and sigma of OrnsteinUhlenbeck

Bridge.b(t,x, P::OrnsteinUhlenbeck) = -P.β*x
Bridge.σ(t, x, P::OrnsteinUhlenbeck) = P.σ
Bridge.a(t, x, P::OrnsteinUhlenbeck) = P.σ^2
Bridge.Γ(t, x, P::OrnsteinUhlenbeck) = inv(P.σ^2)
Bridge.transitionprob(s, x, t, P::OrnsteinUhlenbeck) = Normal(x*exp(-P.β*(t-s)), sqrt((0.5P.σ^2/P.β) *(1-exp(-2*P.β*(t-s)))))
Bridge.constdiff(::OrnsteinUhlenbeck) = true

if ! @isdefined(VOrnsteinUhlenbeck)
struct VOrnsteinUhlenbeck{d}  <: ContinuousTimeProcess{SVector{d,Float64}}
    β
    σ
end
end

# define drift and sigma of VOrnsteinUhlenbeck

Bridge.b(t, x, P::VOrnsteinUhlenbeck) = -P.β*x
Bridge.σ(t, x, P::VOrnsteinUhlenbeck) = P.σ*I
Bridge.bderiv(t, x, P::VOrnsteinUhlenbeck) = -P.β*I
Bridge.σderiv(t, x, P::VOrnsteinUhlenbeck) = 0.0*I
Bridge.a(t, x, P::VOrnsteinUhlenbeck) = P.σ*P.σ'*I
Bridge.constdiff(::VOrnsteinUhlenbeck) = true

kernel(x, a=0.001) = 1/sqrt(2pi*a)* exp(-abs2(x)/(2a))

n = 500
tt = 1.0:1/n:2.0
m = 1000
P = VOrnsteinUhlenbeck{2}(2.0, 1.0)
P1 = OrnsteinUhlenbeck(2.0, 1.0)

u = @SVector [0.0, 0.0]
v = @SVector [0.5, 0.0]
L = @SMatrix [1.0 0.0]
S = @SMatrix [1.0]

convolution(N1::Normal,N2::Normal) = Normal(mean(N1)+mean(N2), sqrt(var(N1)+var(N2)))
EXgivenXpY(X, Y, z) = (z) * var(X) / (var(X) + var(Y))

X = [solve(EulerMaruyama(), u, sample(tt, Wiener{SVector{2,Float64}}()),P).yy[end][1] for i in 1:m]
Xstat = mean(X),var(X)
p1(s, x, t, P::VOrnsteinUhlenbeck) = Normal(x*exp(-P.β*(t-s)), sqrt((0.5P.σ^2/P.β) *(1-exp(-2*P.β*(t-s)))))
PX1 = p1(tt[1],u[1], tt[end], P)
Peta = Normal(0,1.0)
Pobs = convolution(PX1, Peta)
p = pdf(PX1,v[1])
pt= pdf(Normal(0.0, 1.0),v[1])
p2 = pdf(Pobs,v[1])
pt2 = pdf(Normal(0.0, sqrt(2.0)),v[1])

X = solve(EulerMaruyama(), 0.0, sample(tt, Wiener{Float64}()),P1)
@test abs(girsanov(X, P1, Wiener{Float64}()) - llikelihood(X, P1) + llikelihood(X,Wiener{Float64}())) < 0.5

y = EXgivenXpY(PX1, Peta, v[1])

Po = FilterProp(P, tt[1], u, tt[end], v, L, S, Bridge.a(tt[end],v,P))
Po3 = PBridgeProp(P, tt[1], u, (tt[end]-tt[1])/2, 1.2v,tt[end], v, L, S, Bridge.a(tt[end],v,P))

Y = Float64[
let
 X = solve(EulerMaruyama(), (@SVector [0.0, 0.0]), sample(tt, Wiener{SVector{2,Float64}}()),Po)
 X.yy[end][1]*exp(llikelihood(X, Po))*pt2/p2
 end

 for i in 1:m]

Z = Float64[
let
 X = solve(EulerMaruyama(), (@SVector [0.0, 0.0]), sample(tt, Wiener{SVector{2,Float64}}()),Po)
 exp(llikelihood(X, Po))*pt2/p2
 end
 for i in 1:m]

@test abs(mean(Y .- y)*sqrt(m)/std(Y)) < percentile
@test abs(mean(Z .- 1)*sqrt(m)/std(Z)) < percentile

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

Cnames = []
C = []

# BridgeProp
push!(Cnames, "BridgeProp")
Random.seed!(5)
n, m = 200, 1000
T = 2.
ss = range(0, stop=T, length=n)
tau(s, T) = s.*(2-s/T)
#tt = tau(ss, T)
tt = ss

u = 0.5
v = 0.1
a = 0.7
P1 = OrnsteinUhlenbeck(0.8, sqrt(a))
#cs = Bridge.CSpline(tt[1], tt[end], u, u*exp(-P1.β*T), -P1.β*u, -exp(-P1.β*T)*P1.β*u)
#cs = Bridge.CSpline(tt[1], tt[end],  -P1.β*u, -exp(-P1.β*T)*P1.β*u)
h = 0.01

cs = Bridge.CSpline(tt[1], tt[end],
    Bridge.b(tt[1], u, P1),
    Bridge.b(tt[end], v, P1),
    (Bridge.b(tt[2], u + Bridge.b(tt[1], u, P1)*(tt[2]-tt[1]), P1)-Bridge.b(tt[1], u, P1))/(tt[2]-tt[1]), # -P1.β*u*(1-exp(-P1.β*dt))/dt
    (Bridge.b(tt[end], v, P1) - Bridge.b(tt[end-1], v - Bridge.b(tt[end], v, P1)*(tt[end]-tt[end-1]), P1))/(tt[end]-tt[end-1])
)

Po = BridgeProp(P1, tt, (u, v), a, cs)
Z = Float64[
    let
    X = solve(EulerMaruyama(), u, sample(tt, Wiener{Float64}()),Po)
    exp(llikelihood(X, Po))
    end
    for i in 1:m]

p = pdf(transitionprob(0.0, u, T, P1), v)

Pt = Bridge.Ptilde(cs, sqrt(a))
pt = exp(lp(0.0, u, T, v, Pt))
@test pt ≈ exp(lptilde(Po))
push!(C, abs(mean(Z*pt/p .- 1)*sqrt(m)/std(Z*pt/p)))

# Scaling for BridgeProp
push!(Cnames, "ScaledBridgeProp")
Z = Float64[
    let
        X = Bridge.ubridge(sample(tt, Wiener{Float64}()), Po)
        exp(llikelihood(X, Po))
    end
    for i in 1:m]

p = pdf(transitionprob(0.0, u, T, P1), v)

Pt = Bridge.ptilde(Po)
pt = exp(lptilde(Po))
push!(C, abs(mean(Z*pt/p .- 1)*sqrt(m)/std(Z*pt/p)))
#error("end")

# GuidedProp
push!(Cnames, "GuidedProp")
β = 0.8
Ptarget = LinPro(-β, 0.0, sqrt(a))
Po = GuidedProp(Ptarget, tt[1], u, tt[end], v, Pt)

z = Float64[
    let
    W = sample(tt, Wiener{Float64}())
    X = Bridge.bridgeold!(copy(W), W, Po)
    llikelihood(X, Po)
    end
    for i in 1:m]

p2 = pdf(transitionprob(0.0, u, T, Ptarget), v)
p = exp(lp(0.0, u, T, v, Ptarget))
pt = exp(lp(0.0, u, T, v, Pt))
@test p == p2
@test pt ≈ exp(lptilde(Po))
push!(C, abs(mean(exp.(z)*pt/p .- 1)*sqrt(m)/std(exp.(z)*pt/p)))

# DHBridgeProp
push!(Cnames, "DHBridgeProp")

Po3 = DHBridgeProp(P1, tt[1], u, tt[end], v)
Z = Float64[
    let
    X = solve(EulerMaruyama(), u, sample(tt, Wiener{Float64}()),Po3)
    exp(llikelihood(X, Po3))
    end
    for i in 1:m]

p = pdf(transitionprob(0.0, u, T, P1), v)
pt = exp(lptilde(Po3))
@test lptilde(Po3) ≈ Bridge.logpdfnormal(v-u, T*a)
push!(C, abs(mean(Z*pt/p .- 1)*sqrt(m)/std(Z*pt/p)))

# PBridgeProp
push!(Cnames, "PBridgeProp")
tm, vm = 0.5, 0.7
si = 1.
L = 1.
Po2 = PBridgeProp(P1, tt[1], u, tm, vm, tt[end], v, L, si^2, a, cs)
Z2 = Float64[
    let
    X = solve(EulerMaruyama(), u, sample(tt, Wiener{Float64}()),Po2)
    exp(llikelihood(X, Po2))
    end
    for i in 1:m]

f(x) = pdf(transitionprob(0.0, u, tm, P1), x)*pdf(transitionprob(tm,x,T, P1), v)*kernel.(x-vm,si^2)
ft(x) = exp(Bridge.lp(0.0, u, tm, x, Pt) + Bridge.lp(tm,x,T, v, Pt))*kernel.(x-vm,si^2)
p2 = sum(map(f,range(-20, stop=20, length=1001)))*40/1000
pt2 = exp(Bridge.lptilde(Po2))
@test pt2 ≈ sum(map(ft,range(-20, stop=20, length=1001)))*40/1000
push!(C, abs(mean(Z2*pt2/p2 .- 1)*sqrt(m)/std(Z2*pt2/p2)))

# GuidedProp
push!(Cnames, "GuidedProp")
Ptarget = Bridge.Ptilde(cs, sqrt(a))
Pt = LinPro(-β, 0.2, sqrt(a))
Po = GuidedProp(Ptarget, tt[1], u, tt[end], v, Pt)

z = Float64[
    let
    W = sample(tt, Wiener{Float64}())
    X = Bridge.bridgeold!(copy(W), W, Po)
    llikelihood(X, Po)
    end
    for i in 1:m]

p2 = pdf(transitionprob(0.0, u, T, Ptarget), v)
p = exp(lp(0.0, u, T, v, Ptarget))
pt = exp(lp(0.0, u, T, v, Pt))
@test p ≈ p2
@test pt ≈ exp(lptilde(Po))
push!(C, abs(mean(exp.(z)*pt/p .- 1)*sqrt(m)/std(exp.(z)*pt/p)))

# GuidedProposal
push!(Cnames, "GuidedProposal")
Ptarget = Bridge.Ptilde(cs, sqrt(a))
Pt = LinPro(-β, 0.2, sqrt(a))
GP = Bridge.GuidedBridge(tt, Ptarget, Pt, v)

@test norm(GP.H♢-[inv(Bridge.H(t, T, Pt)) for t in tt], Inf) < 1e-5
@test norm(GP.V-[Bridge.V(t, T, v, Pt) for t in tt], Inf) < 1e-5

z = Float64[
    let
    W = sample(tt, Wiener{Float64}())
    X = copy(W)
    Bridge.solve!(Euler(), X, u, W, GP)
    llikelihood(LeftRule(), X, GP)
    end
    for i in 1:m]

p2 = pdf(transitionprob(0.0, u, T, Ptarget), v)
p = exp(lp(0.0, u, T, v, Ptarget))
pt = exp(lp(0.0, u, T, v, Pt))
@test p ≈ p2
@test pt ≈ exp(lptilde(Po))
push!(C, abs(mean(exp.(z)*pt/p .- 1)*sqrt(m)/std(exp.(z)*pt/p)))

println(Cnames)
println(C)

@test  all(C .< percentile)


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# Contributed by Adam Beckmeyer. Based on code by Jarret Revels, Alex
# Arslan, Michal Stransky, Jens Adam.

using Distributed

@everywhere struct Node
    l::Union{Node,Nothing}
    r::Union{Node,Nothing}
end

@everywhere make(n) =
    n === 0 ? Node(nothing, nothing) : Node(make(n-1), make(n-1))

@everywhere check(node) =
    node.l === nothing ? 1 : 1 + check(node.l) + check(node.r)

function binary_trees(io, n)
    write(io, "stretch tree of depth $(n+1)\t check: $(check(make(n+1)))\n")

    long_tree = make(n)

    d = 4
    while d <= n
        niter = 1 << (n - d + 4)
        c = @distributed (+) for _ in 1:niter
            check(make(d))
        end#for
        write(io, "$niter\t trees of depth $d\t check: $c\n")
        d += 2
    end

    write(io, "long lived tree of depth $n\t check: $(check(long_tree))\n")
end#function

isinteractive() || binary_trees(stdout, parse(Int, ARGS[1]))

max_stars_repo_name: bpecsek/Programming-Language-Benchmarks

max_issues_repo_name: bpecsek/Programming-Language-Benchmarks

max_forks_repo_name: bpecsek/Programming-Language-Benchmarks


text:
using Test
using Primes
using FFTW

import OpenCL.cl

import CLFFT
const clfft = CLFFT

#macro throws_pred(ex) FactCheck.throws_pred(ex) end

const TOLERANCE = 1e-3

function allclose(x::AbstractArray{T}, y::AbstractArray{T}; rtol=1e-5, atol=1e-8) where {T}
    @assert length(x) == length(y)
    @inbounds begin
        for i in length(x)
            xx, yy = x[i], y[i]
            if !(isapprox(xx, yy; rtol=rtol, atol=atol))
                return false
            end
        end
    end
    return true
end

#Note: these functions have been adapted from the clFFT tests (buffer.h)
function floats_about_equal(a::T, b::T) where {T<:clfft.clfftNumber}
    tol  = convert(T, TOLERANCE)
    comp = 1e-5
    if T <: Real
        if abs(a) < comp && abs(b) < comp
            return true
        end
        return abs(a - b) > abs(a * tol) ? false : true
    else
        return (floats_about_equal(real(a), real(b)) &&
                floats_about_equal(imag(b), imag(b)))
    end
end

function allclose_clfft(x::AbstractArray{T}, y::AbstractArray{T}) where {T<:clfft.clfftNumber}
    @assert length(x) == length(y)
    @inbounds begin
        for i in length(x)
            xx, yy = x[i], y[i]
            if !(floats_about_equal(xx, yy))
                return false
            end
        end
    end
    return true
end

@testset "CLFFT.jl" begin
    @testset "Version" begin
        @test isa(CLFFT.version(), VersionNumber)
        v = clfft.version()
        @test v.major >= 2
        @test v.minor >= 1
        @test v.patch >= 0
    end

    @testset "Example FFT Single" begin
        for N in (2^8,)# 3^7, 5^6]
            X = rand(ComplexF32, N)
            fftw_X = fft(X)
            for device in cl.devices()
                @testset "$(device[:name])" begin
                    ctx = cl.Context(device)
                    queue = cl.CmdQueue(ctx)
                    bufX = cl.Buffer(ComplexF32, ctx, :copy, hostbuf=X)
                    p = clfft.Plan(ComplexF32, ctx, size(X))
                    clfft.set_layout!(p, :interleaved, :interleaved)
                    clfft.set_result!(p, :inplace)
                    clfft.bake!(p, queue)

                    @test clfft.context(p) == ctx
                    @test clfft.precision(p) == :single
                    @test clfft.layout(p) == (:interleaved, :interleaved)
                    @test clfft.result(p) == :inplace
                    @test clfft.dim(p) == 1
                    @test length(clfft.lengths(p)) == 1
                    @test clfft.lengths(p)[1] == length(X)
                    @test clfft.transpose_result(p) == false

                    @test clfft.scaling_factor(p, :forward) == Float32(1.0)
                    @test clfft.scaling_factor(p, :backward) == Float32(1.0 / length(X))
                    @test clfft.batchsize(p) == 1

                    clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
                    # read is blocking (waits on pending event for result)
                    R = cl.read(queue, bufX)
                    @test allclose(R, fftw_X; rtol=1e-2, atol=1e-3)
                    @test allclose_clfft(R, fftw_X)
                end
            end
        end
    end

    @testset "Example FFT Double" begin
        for N in (2^7,)# 3^6, 5^5]
            X = rand(ComplexF64, N)
            fftw_X = fft(X)
            for device in cl.devices()
                @testset "$(device[:name])" begin
                    ctx = cl.Context(device)
                    queue = cl.CmdQueue(ctx)
                    bufX = cl.Buffer(ComplexF64, ctx, :copy, hostbuf=X)
                    p = clfft.Plan(ComplexF64, ctx, size(X))
                    clfft.set_layout!(p, :interleaved, :interleaved)
                    clfft.set_result!(p, :inplace)
                    clfft.bake!(p, queue)

                    @test clfft.context(p) == ctx
                    @test clfft.precision(p) == :double
                    @test clfft.layout(p) == (:interleaved, :interleaved)
                    @test clfft.result(p) == :inplace
                    @test clfft.dim(p) == 1
                    @test length(clfft.lengths(p)) == 1
                    @test clfft.lengths(p)[1] == length(X)
                    @test clfft.transpose_result(p) == false

                    @test clfft.scaling_factor(p, :forward) == Float32(1.0)
                    @test clfft.scaling_factor(p, :backward) == Float32(1.0 / length(X))
                    @test clfft.batchsize(p) == 1

                    clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
                    R = cl.read(queue, bufX)
                    @test allclose(R, fftw_X; rtol=1e-2, atol=1e-3)
                    @test allclose_clfft(R, fftw_X)
                end
            end
        end
    end

    @testset "2D FFT Inplace" begin
        transform_sizes = (2^6,)#3^4, 5^3
        for N in transform_sizes
            for M in transform_sizes
                X = rand(ComplexF32, (N, M))
                fftw_X = fft(X)
                for device in cl.devices()
                    @testset "$(device[:name])" begin
                        ctx = cl.Context(device)
                        queue = cl.CmdQueue(ctx)
                        p = clfft.Plan(ComplexF32, ctx, X)
                        clfft.bake!(p, queue)
                        bufX = cl.Buffer(ComplexF32, ctx, :copy, hostbuf=X)
                        clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
                        R = reshape(cl.read(queue, bufX), size(X))
                        @test allclose(R, fftw_X; rtol=1e-2, atol=1e-3)
                        @test allclose_clfft(R, fftw_X)
                        GC.gc()
                    end
                end
            end
        end
    end

    @testset "3D FFT Inplace" begin
        N = 64
        X = rand(ComplexF32, (N, N, N))
        fftw_X = fft(X)
        for device in cl.devices()
            @testset "$(device[:name])" begin
                ctx = cl.Context(device)
                queue = cl.CmdQueue(ctx)
                p = clfft.Plan(ComplexF32, ctx, X)
                clfft.bake!(p, queue)
                bufX = cl.Buffer(ComplexF32, ctx, :copy, hostbuf=X)
                clfft.enqueue_transform(p, :forward, [queue], bufX, nothing)
                R = reshape(cl.read(queue, bufX), size(X))
                @test allclose(R, fftw_X; rtol=1e-2, atol=1e-3)
                @test allclose_clfft(R, fftw_X)
                GC.gc()
            end
        end
    end
end


text:
"""
    resonance_eigenproblem(simulation, k, ka=0, kb=0) -> A,B,σ

generate matrices from `simulation` to solve {∇²+εk²}ψ=0 in the form `A*v=λB*b` with shift `σ`
"""
function resonance_eigenproblem(sim::Simulation, k::Number, ka::Number=0, kb::Number=0)
    L = Laplacian(sim,k,ka,kb)
    return L.D, -L.M, k^2
end

"""
    resonance_eigenproblem(sim, k, ka=0, kb=0; η=0) -> A,B,σ

generate matrices from `simulation` to solve {∇²+εk²}u=-ηFk²u in the form `A*v=λB*b` with shift `σ
"""
function cf_eigenproblem(sim::Simulation, k::Number, ka::Number=0, kb::Number=0; η::Number=0)
    L = Laplacian(sim,k,ka,kb)
    return L.D+L.M*k^2, -L.Mf*k^2, η
end


text:
using Test, LinearAlgebra, BasisFunctions, FrameFun

dict1 = Fourier(10)
dict2 = Fourier(10)^2
plat1 = platform(dict1)
plat2 = platform(dict2)
par1 = 10
par2 = (10,10)

all_dicts = (dict1, dict2)

@testset "samplingparameter" begin
    # adaptive approximations has no samplingparameter
    for (a,sol) in zip(all_dicts ,(10,(10,10),10,(10,10)))
        @test samplingparameter(a)==sol
    end
    @test samplingparameter(plat1, par1) == 10
    @test samplingparameter(plat2, par2) == (10,10)
end

@testset "interpolation_grid" begin
    @test interpolation_grid(plat1, par1) == interpolation_grid(dict1)
    @test interpolation_grid(plat2, par2) == interpolation_grid(dict2)

    ap = approximationproblem(dict1)
    @test !FrameFun.ap_hasproperty(ap, interpolation_grid)
    @test interpolation_grid(ap) == interpolation_grid(dict1)
    @test FrameFun.ap_hasproperty(ap, interpolation_grid)
    @test FrameFun.cache(ap, interpolation_grid) == interpolation_grid(dict1)
end

@testset "oversampling_grid" begin
    @test oversampling_grid(dict1) == FourierGrid(10)
    @test oversampling_grid(dict1, 11) == FourierGrid(11)
    @test oversampling_grid(dict2) == FourierGrid(10)^2
    @test oversampling_grid(dict2, (11,11)) == FourierGrid(11)^2
    @test oversampling_grid(plat1, par1) == FourierGrid(10)

    ap = FrameFun.approximationproblem(dict1, 15)
    @test !FrameFun.ap_hasproperty(ap, oversampling_grid)
    @test oversampling_grid(ap) == FourierGrid(15)
    @test FrameFun.ap_hasproperty(ap, oversampling_grid)
    @test FrameFun.cache(ap, oversampling_grid) == FourierGrid(15)
end

@testset "sampling_grid" begin
    @test sampling_grid(dict1) == FourierGrid(10)
    @test sampling_grid(dict2) == FourierGrid(10)^2
    @test sampling_grid(plat1, par1) == FourierGrid(10)
    @test sampling_grid(plat2, par2) == FourierGrid(10)^2
end

@testset "platform_grid" begin
    @test platform_grid(dict1; samplingstyle=GridStyle(), grid=FourierGrid(1)) == FourierGrid(1)
    @test platform_grid(dict2; samplingstyle=GridStyle(), grid=FourierGrid(1)^2) == FourierGrid(1)^2
    @test platform_grid(plat1, par1, grid=FourierGrid(1)) == FourierGrid(1)
    @test sampling_grid(dict2; samplingstyle=GridStyle(),grid=FourierGrid(1)^2)==FourierGrid(1)^2
end

@testset "sampling_operator" begin
    op1 = sampling_operator(dict1)
    op2 = sampling_operator(dict2)
    op3 = sampling_operator(plat1, par1)
    op4 = sampling_operator(plat2, par2)

    @test grid(dest(op1)) == FourierGrid(10)
    @test grid(dest(op3)) == FourierGrid(10)

    f2 = (x,y)->exp(x*y)
    @test tensorproduct(op1,op1)*f2≈op2*f2
    @test tensorproduct(op3,op3)*f2≈op4*f2
end

@testset "discretemeasure" begin
    @test discretemeasure(dict1) == discretemeasure(sampling_grid(dict1))
    @test discretemeasure(dict2)== discretemeasure(sampling_grid(dict2))
    @test discretemeasure(plat1, par1)== discretemeasure(sampling_grid(plat1, par1))
    @test discretemeasure(plat2, par2)== discretemeasure(sampling_grid(plat2, par2))
end

@testset "measure" begin
    @test measure(platform(dict1)) == measure(dict1)
    @test measure(platform(dict2)) == measure(dict2)
end

@testset "azdual" begin
    @test operator(azdual(dict1))≈operator(azdual(plat1, par1))
    op = azdual(dict2)
    @test op isa TensorProductDict
    @test ≈(operator.(components(op))...)
    op = azdual(plat2, par2)
    @test op isa TensorProductDict
    @test ≈(operator.(components(op))...)
    @test operator(component(azdual(dict2),1))≈operator(azdual(dict1))

    @test operator(azdual(dict1;samplingstyle=GramStyle())) ≈
        operator(component(azdual(dict2;samplingstyle=ProductSamplingStyle(GramStyle(),GramStyle())),1))
    @test operator(azdual(dict1;samplingstyle=GramStyle())) ≈
        operator(component(azdual(plat2, par2;samplingstyle=ProductSamplingStyle(GramStyle(),GramStyle())),1))
    @test operator(azdual(dict1;samplingstyle=GramStyle())) ≈
        operator(component(azdual(dict2;samplingstyle=ProductSamplingStyle(GramStyle(),GramStyle())),1))
end

@testset "discretization" begin
    op1 = discretization(plat1, par1)
    op2 = discretization(plat2, par2)
    op3 = discretization(dict1)
    op4 = discretization(dict2)

    @test op2 isa TensorProductOperator
    @test op4 isa TensorProductOperator

    @test op1 ≈ op3
    @test op2 ≈ op4

    f1 = exp
    f2 = (x,y)->exp(x*y)
    op3, b3 = full_discretization(f1, plat1, par1)
    op4, b4 = full_discretization(f2, plat2, par2)
    op5, b5 = full_discretization(f1, dict1)
    op6, b6 = full_discretization(f2, dict2)

    @test op4 isa TensorProductOperator
    @test op6 isa TensorProductOperator

    @test op3 ≈ op5
    @test op4 ≈ op6

    @test b3≈b5
    @test b4≈b6
end

@testset "dualdiscretization" begin
    op1 = dualdiscretization(plat1, par1)
    op2 = dualdiscretization(plat2, par2)
    op3 = dualdiscretization(dict1)
    op4 = dualdiscretization(dict2)

    @test op2 isa TensorProductOperator
    @test op4 isa TensorProductOperator

    @test op1 ≈ op3
    @test op2 ≈ op4
end

@testset "solver" begin
    op1 = solver(dict1)
    op2 = solver(dict2)
    op3 = solver(plat1, par1)
    op4 = solver(plat2, par2)

    @test op2 isa TensorProductOperator
    @test op4 isa TensorProductOperator

    @test op1 ≈ op3
    @test op2 ≈ op4

    @test component(op2,1) ≈ op1 ≈ component(op2,2)

    for STYLE in (InverseStyle, DirectStyle, DualStyle, IterativeStyle,AZStyle,AZSmoothStyle)
        op1 = solver(dict1;solverstyle=STYLE())
        op2 = solver(dict2;solverstyle=ProductSolverStyle(STYLE(),STYLE()))
        op3 = solver(plat1, par1;solverstyle=STYLE())
        op4 = solver(plat2, par2;solverstyle=ProductSolverStyle(STYLE(),STYLE()))

        @test op2 isa TensorProductOperator
        @test op4 isa TensorProductOperator

        @test op1 ≈ op3
        @test op2 ≈ op4

        @test component(op2,1) ≈ op1 ≈ component(op2,2)
    end
end

@testset "AZ_A, AZ_Z, AZ_Zt, plungeoperator, firstAZstepoperator, plungerank" begin
    A1 = AZ_A(dict1)
    A2 = AZ_A(dict2)
    @test tensorproduct(A1,A1) ≈ A2

    Z1 = AZ_Z(dict1)
    Z2 = AZ_Z(dict2)
    @test tensorproduct(Z1,Z1) ≈ Z2

    Zt1 = AZ_Zt(dict1)
    Zt2 = AZ_Zt(dict2)
    @test tensorproduct(Zt1,Zt1) ≈ Zt2

    @test Zt2 ≈ Z2'

    P1 = plungeoperator(dict1)
    P2 = plungeoperator(dict2)
    @test 1+norm(tensorproduct(P1,P1) -P2)≈1
    @test I-A2*Zt2 ≈ P2

    M1 = firstAZstepoperator(dict1)
    M2 = firstAZstepoperator(dict2)
    @test 1+norm(tensorproduct(M1,M1) - M2)≈1
    @test 1+norm(A2*P2 - M2)≈1

    @test plungerank(dict1) <= 1
    @test plungerank(dict2) <= 1
end

@testset begin "ParameterPath"
    p = IncrementalCartesianParameterPath{2}()
    P = ExtensionFramePlatform(WeightedSumPlatform(platform(ChebyshevT(10)^2), (x,y)->1.,
                (x,y)->sqrt(x^2+y^2)), .9 .* UnitDisk())
    path = HierarchyPath(p,ProductPath(p,p))
    paramP = parametrizedplatform(P, path)

    paramPparam=45
    Pparam = path[paramPparam]

    @test SamplingStyle(P)==SamplingStyle(paramP)
    @test SolverStyle(P, SamplingStyle(P)) == SolverStyle(paramP, SamplingStyle(P))
    @test sampling_grid(P,Pparam)≈sampling_grid(paramP, paramPparam)
    @test dimensions(dictionary(P,Pparam)) == dimensions(dictionary(paramP,paramPparam))
    @test dimensions(azdual(P,Pparam))==dimensions(azdual(paramP, paramPparam))
    @test AZ_A(P,Pparam)≈AZ_A(paramP, paramPparam)
    @test AZ_Zt(P,Pparam)≈AZ_Zt(paramP, paramPparam)
    @test AZ_A(P,Pparam)≈AZ_A(paramP, paramPparam)
    @test typeof(solver(P,Pparam))==typeof(solver(paramP,paramPparam))
end


text:
function adjoint_Jacobi_field(
    M::NONMUTATINGMANIFOLDS,
    p::Real,
    q::Real,
    t::Real,
    X::Real,
    β=βdifferential_geodesic_startpoint,
)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M, p, q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (inner.(Ref(M), Ref(x), Ref(X), Θ)) .*
        (β.(B.data.eigenvalues, Ref(t), Ref(distance(M, p, q)))) .* V,
    )[1]
    return Y
end
function jacobi_field(
    M::NONMUTATINGMANIFOLDS,
    p::Real,
    q::Real,
    t::Real,
    X::Real,
    β=βdifferential_geodesic_startpoint,
)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M, p, q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    Y = zero_tangent_vector(M, p)
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (inner.(Ref(M), Ref(p), Ref(X), V)) .*
        (β.(B.data.eigenvalues, Ref(t), Ref(distance(M, p, q)))) .* Θ,
    )[1]
    return Y
end
function grad_TV2(M::NONMUTATINGMANIFOLDS, q, p::Int=1)
    c = mid_point(M, q[1], q[3], q[2]) # nearest mid point of x and z to y
    d = distance(M, q[2], c)
    innerLog = -log(M, c, q[2])
    X = [zero_tangent_vector(M, q[i]) for i in 1:3]
    if p == 2
        X[1] = adjoint_differential_geodesic_startpoint(M, q[1], q[3], 1 / 2, innerLog)
        X[2] = -log(M, q[2], c)
        X[3] = adjoint_differential_geodesic_endpoint(M, q[1], q[3], 1 / 2, innerLog)
    else
        if d > 0 # gradient case (subdifferential contains zero, see above)
            X[1] = adjoint_differential_geodesic_startpoint(
                M, q[1], q[3], 1 / 2, innerLog / (d^(2 - p))
            )
            X[2] = -log(M, q[2], c) / (d^(2 - p))
            X[3] = adjoint_differential_geodesic_endpoint(
                M, q[1], q[3], 1 / 2, innerLog / (d^(2 - p))
            )
        end
    end
    return X
end
function prox_TV2(::NONMUTATINGMANIFOLDS, λ, pointTuple::Tuple{T,T,T}, p::Int=1) where {T}
    w = @SVector [1.0, -2.0, 1.0]
    x = SVector(pointTuple)
    if p == 1 # Theorem 3.5 in Bergmann, Laus, Steidl, Weinmann, 2014.
        sr_dot_xw = sym_rem(sum(x .* w))
        m = min(λ, abs(sr_dot_xw) / (dot(w, w)))
        s = sign(sr_dot_xw)
        return sym_rem.(x .- m .* s .* w)
    elseif p == 2 # Theorem 3.6 ibd.
        t = λ * sym_rem(sum(x .* w)) / (1 + λ * dot(w, w))
        return sym_rem.(x - t .* w)
    else
        throw(
            ErrorException(
                "Proximal Map of TV2(Circle,λ,pT,p) not implemented for p=$(p) (requires p=1 or 2)",
            ),
        )
    end
end
function prox_TV2(M::PowerManifold{𝔽,N}, λ, x, p::Int=1) where {𝔽,N<:NONMUTATINGMANIFOLDS}
    power_size = power_dimensions(M)
    R = CartesianIndices(power_size)
    d = length(size(x))
    minInd = first(R).I
    maxInd = last(R).I
    y = deepcopy(x)
    for k in 1:d # for all directions
        ek = CartesianIndex(ntuple(i -> (i == k) ? 1 : 0, d)) #k th unit vector
        for l in 0:2
            for i in R # iterate over all pixel
                if (i[k] % 3) == l
                    JForward = i.I .+ ek.I #i + e_k
                    JBackward = i.I .- ek.I # i - e_k
                    if all(JForward .<= maxInd) && all(JBackward .>= minInd)
                        (y[M, JBackward...], y[M, i.I...], y[M, JForward...]) = prox_TV2(
                            M.manifold,
                            λ,
                            (y[M, JBackward...], y[M, i.I...], y[M, JForward...]),
                            p,
                        )
                    end
                end # if mod 3
            end # i in R
        end # for mod 3
    end # directions
    return y
end


text:
module SubDiff
using ValidatedNumerics;
include("sti.jl")

function productLoSubDifferential{T}(aInterval :: Interval{T}, xLo :: T, xHi :: T)
    maxOfFactor = dMaxHiLo(aInterval, sti.pos(aInterval.hi)*sti.pos(xLo), sti.neg(aInterval.lo)*sti.pos(xHi))
    result = sti.pos(aInterval.lo) * dXLoNeg(xLo) + sti.neg(aInterval.hi) * dXHiNeg(xHi) - maxOfFactor
    result
end

function productHiSubDifferential{T}(aInterval :: Interval{T}, xLo :: T, xHi :: T)
    maxOfFactor = dMaxHiHi(aInterval, sti.pos(aInterval.hi)*sti.pos(xHi), sti.neg(aInterval.lo)*sti.pos(xLo))
    result = - sti.pos(aInterval.lo) * dXHiNeg(xHi) - sti.neg(aInterval.hi) * dXLoNeg(xLo) + maxOfFactor
    result
end

function dXLoPos(x)
    if x < 0
        @interval(0, 0)
    elseif x == 0
        @interval(-0.5, 0)
    else
        @interval(-1, 0)
    end
end

function dXLoNeg(x)
    if x < 0
        @interval(0, 1)
    elseif x == 0
        @interval(0, 0.5)
    else
        @interval(0)
    end
end

function dXHiPos(x)
    if x < 0
        @interval(0)
    elseif x == 0
        @interval(0, 0.5)
    else
        @interval(0, 1)
    end
end

function dXHiNeg(x)
    if x < 0
        @interval(-1, 0)
    elseif x == 0
        @interval(-0.5, 0)
    else
        @interval(0)
    end
end

function dMaxHiLo(c, left, right)
    if left > right
        @interval(-sti.pos(c.hi), 0)
    elseif left == right
        0.5*@interval(-sti.pos(c.hi), sti.neg(c.lo))
    else
        @interval(0, sti.neg(c.lo))
    end
end

function dMaxHiHi(c, left, right)
    if left > right
        @interval(0, sti.pos(c.hi))
    elseif left == right
        0.5*@interval(-sti.neg(c.lo), sti.pos(c.hi))
    else
        @interval(-sti.neg(c.lo), 0)
    end
end

end


text:
#
# This is remote terminal application. It's used for remote
# configuring and statistics collecting of server. It supports
# all commands, which are declared in RpcApi.jl module.
#
# @author DeadbraiN
#
include("global/ImportFolders.jl")

import Client
import Creature
import Connection
import Config
importall Terminal
#
# Answer function. It's called every time on request from server
# @param ans Server answer
#
function onAnswer(ans::Connection.Answer)
  global answer = ans.data
  sleep(0.3)
  println(answer)
end
#
# Application entry point function. Adds listeners for
# server responses and initializes Terinal module.
#
function main()
  local cfg::Config.ConfigData = Config.create()
  # TODO: ip should be get from config
  local con = Client.create(ip"127.0.0.1", cfg.conServerPort)

  if !Client.isOk(con) || Config.isEmpty(cfg) quit() end
  Event.on(con.observer, Client.EVENT_AFTER_REQUEST, onAnswer)
  Terminal.init(con)
end
#
# App start point
#
main()

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"This script will run the perturbation function on each case of increasing size
(determined by the parameter num_cases) and attempt to solve. If a solution times
out (at 3000 iterations) the result is written to the unsolved/ folder. Otherwise,
the output of each run is fed to result_dict, which contains tuples indicating the
(result_pert_loss, result_pert_cost)."

using Ipopt
using PowerModelsPrivacyPreserving
using PowerModels
using Distributions
using Random
using JuMP

const PMs = PowerModels
const PMPP = PowerModelsPrivacyPreserving

"https://stackoverflow.com/questions/48195775/how-to-pretty-print-nested-dicts-in-julia
Adapted by EV to write to file"
function pretty_print_to_file(io, d::Dict, pre=1)
    for (k,v) in d
        if typeof(v) <: Dict
            s = "$(repr(k)) => "
            write(io, join(fill(" ", pre)) * s * "\n")
            pretty_print_to_file(io, v, pre+1+length(s))
        else
            write(io, join(fill(" ", pre)) * "$(repr(k)) => $(repr(v))" * "\n")
        end
    end
    nothing
end

function check_dataset_perturbation(test_directory, output_directory, filename, α, β, ϵ, λ)
    # ipopt = Ipopt.Optimizer
    optimizer = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "max_cpu_time" => 600.0)
    data_unpert = parse_file(string(test_directory, filename))
    data_min_loss = deepcopy(data_unpert)
    data_min_cost = deepcopy(data_unpert)

    "this is the canonical OPF problem with generation cost minimization"
    result_unpert_cost = PMPP.run_ac_opf_cost(data_unpert, optimizer)
    PMPP.calculate_losses!(result_unpert_cost, data_unpert)
    "store faithfulness info"
    data_min_cost["cost"] = Dict()
    data_min_cost["cost"]["value"] = result_unpert_cost["objective"]
    data_min_cost["cost"]["beta"] = β

    "this variant of the OPF problem minimizes grid losses instead of generation cost"
    result_unpert_loss = PMPP.run_ac_opf_loss(data_unpert, optimizer)
    PMPP.calculate_losses!(result_unpert_loss, data_unpert)
    "store faithfulness info"
    data_min_loss["loss"] = Dict()
    data_min_loss["loss"]["value"] = result_unpert_loss["totalloss"]
    data_min_loss["loss"]["beta"] = β

    # Add impedance perturbation to both data dictionaries.
    data_pert_min_loss = PMPP.create_impedance_perturbation(data_min_loss, α, ϵ, λ)
    data_pert_min_cost = PMPP.create_impedance_perturbation(data_min_cost, α, ϵ, λ)

    # 1) Run solver for perturbed loss
    result_pert_loss = PMPP.run_opf_variable_impedance_loss(data_pert_min_loss, optimizer)
    PMPP.calculate_losses!(result_pert_loss, data_pert_min_loss)
    PMPP.overwrite_impedances_in_data!(result_pert_loss, data_pert_min_loss)
    # @assert result_pert_loss["termination_status"] == PMs.LOCALLY_SOLVED

    # Check how the r and x parameters have been perturbed and postprocessed
    output_r_values_loss = Dict()
    for (l, branch) in result_pert_loss["solution"]["branch"]
        output_r_values_loss[l] = Dict()

        output_r_values_loss[l]["r_original"] = data_unpert["branch"][l]["br_r"]
        output_r_values_loss[l]["r_pert"] = real(1 / (branch["g"] + branch["b"]im))
        output_r_values_loss[l]["r_pert_ratio"] = output_r_values_loss[l]["r_pert"] / output_r_values_loss[l]["r_original"]

        output_r_values_loss[l]["x_original"] = data_unpert["branch"][l]["br_x"]
        output_r_values_loss[l]["x_pert"] = imag(1 / (branch["g"] + branch["b"]im))
        output_r_values_loss[l]["x_pert_ratio"] = output_r_values_loss[l]["x_pert"] / output_r_values_loss[l]["x_original"]

        output_r_values_loss[l]["x_shunt_original"] = -1 / data_unpert["branch"][l]["b_to"]
        output_r_values_loss[l]["x_shunt_pert"] = -1 / branch["b_shunt"]  # Ignore g when calculating x_shunt
        output_r_values_loss[l]["x_shunt_pert_ratio"] = output_r_values_loss[l]["x_shunt_pert"] / output_r_values_loss[l]["x_shunt_original"]
    end

    # Handle writing results to file based on success criteria
    if (result_pert_loss["termination_status"] == PMs.LOCALLY_SOLVED)
        result_directory = "pert_min_loss/"
    elseif (result_pert_loss["termination_status"] == PMs.TIME_LIMIT)
        result_directory = "timed_out_min_loss/"
    else
        result_directory = "failed_min_loss/"
    end
    open(output_directory * result_directory * filename[1:length(filename) - 2] * "_result.txt", "w") do io
        pretty_print_to_file(io, result_pert_loss)
    end
    open(output_directory * result_directory * filename[1:length(filename) - 2] * "_data.m", "w") do io
        PMs.export_matpower(io, data_pert_min_loss)
    end

    # Write r_pert to file
    open(output_directory * result_directory * filename[1:length(filename) - 2] * "_rx_ratios.txt", "w") do io
        pretty_print_to_file(io, output_r_values_loss)
    end

    # 2) Run solver for perturbed cost
    result_pert_cost = PMPP.run_opf_variable_impedance_cost(data_pert_min_cost, optimizer)
    PMPP.calculate_losses!(result_pert_cost, data_pert_min_cost)
    PMPP.overwrite_impedances_in_data!(result_pert_cost, data_pert_min_cost)
    # @assert result_pert_cost["termination_status"] == PMs.LOCALLY_SOLVED

    # Check how the r and x parameters have been perturbed and postprocessed
    output_r_values_cost = Dict()
    for (l, branch) in result_pert_cost["solution"]["branch"]
        output_r_values_cost[l] = Dict()

        output_r_values_cost[l]["r_original"] = data_unpert["branch"][l]["br_r"]
        output_r_values_cost[l]["r_pert"] = real(1 / (branch["g"] + branch["b"]im))
        output_r_values_cost[l]["r_pert_ratio"] = output_r_values_cost[l]["r_pert"] / output_r_values_cost[l]["r_original"]

        output_r_values_cost[l]["x_original"] = data_unpert["branch"][l]["br_x"]
        output_r_values_cost[l]["x_pert"] = imag(1 / (branch["g"] + branch["b"]im))
        output_r_values_cost[l]["x_pert_ratio"] = output_r_values_cost[l]["x_pert"] / output_r_values_cost[l]["x_original"]

        output_r_values_cost[l]["x_shunt_original"] = -1 / data_unpert["branch"][l]["b_to"]
        output_r_values_cost[l]["x_shunt_pert"] = -1 / branch["b_shunt"]  # Ignore g when calculating x_shunt
        output_r_values_cost[l]["x_shunt_pert_ratio"] = output_r_values_cost[l]["x_shunt_pert"] / output_r_values_cost[l]["x_shunt_original"]
    end

    # Handle writing results to file based on success criteria
    if (result_pert_cost["termination_status"] == PMs.LOCALLY_SOLVED)
        result_directory = "pert_min_cost/"
    elseif (result_pert_cost["termination_status"] == PMs.TIME_LIMIT)
        result_directory = "timed_out_min_cost/"
    else
        result_directory = "failed_min_cost/"
    end
    open(output_directory * result_directory * filename[1:length(filename) - 2] * "_result.txt", "w") do io
        pretty_print_to_file(io, result_pert_loss)
    end
    open(output_directory * result_directory * filename[1:length(filename) - 2] * "_data.m", "w") do io
        PMs.export_matpower(io, data_pert_min_loss)
    end

    # Write r_pert to file
    open(output_directory * result_directory * filename[1:length(filename) - 2] * "_rx_ratios.txt", "w") do io
        pretty_print_to_file(io, output_r_values_cost)
    end
end

"Set the variable num_cases to determine how many cases to solve"
num_cases = 37
start_case = 1
start_index = 1

# Make all directories for outputs
test_directory = "test/data/pglib_tests/"
output_directory = "examples/test_perturbation_outputs/"
try
    mkdir(output_directory)
catch y
    println("Output folder already exists, continuing")
end

α = 0.01
β = 0.5
ϵ = 1
λ = 30
for run_index = start_index:100
    run_output_directory = output_directory * string(run_index) * "/"
    try
        mkdir(run_output_directory)
    catch y
        println("Iteration folder already exists, continuing")
    end

    try
        mkdir(string(run_output_directory, "pert_min_loss"))
        mkdir(string(run_output_directory, "pert_min_cost"))
        mkdir(string(run_output_directory, "failed_min_loss"))
        mkdir(string(run_output_directory, "failed_min_cost"))
        mkdir(string(run_output_directory, "timed_out_min_loss"))
        mkdir(string(run_output_directory, "timed_out_min_cost"))
    catch y
        println("Output subfolders already exist, continuing")
    end

    # Sort the list of cases by size
    sorted_directory = sort(
        readdir(test_directory),
        by = f -> parse(Int, strip(split(f, "_")[3][5:end], ['w', 'o', 'p', 's']))
    )
    for filename in sorted_directory[start_case: num_cases]
        println("Testing ", filename)
        check_dataset_perturbation(test_directory, run_output_directory, filename, α, β, ϵ, λ)
    end
    global start_case = 1
end


text:
# define global
x = 0

# define local
let
    local x = 1
    let
        local x = 2
    end
    println(x) # local x
end
println(x) # global x

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
using Setfield

const cellml_ns(xml::EzXML.Document) = namespace(root(xml))
const cellml_ns(node::EzXML.Node) = namespace(node)
const mathml_ns = "http://www.w3.org/1998/Math/MathML"

# create_var(x) = Num(Sym{Real}(Symbol(x))).val  # Num(Variable(Symbol(x))).val
# create_var(x, iv) = Num(Sym{Symbolics.FnType{Tuple{Any},Real}}(Symbol(x)))(iv).val  # Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(Symbol(x)))(iv).val

create_var(x) = Num(Variable(Symbol(x))).val
create_var(x, iv) = Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(Symbol(x)))(iv).val

function create_param(x)
  p = Sym{Real}(Symbol(x))
  ModelingToolkit.toparam(p)
end

# function create_var(x)
#     @variables v
#     ModelingToolkit.rename(v.val, Symbol(x))
# end
#
# function create_var(x, iv)
#     @variables v(iv)
#     ModelingToolkit.rename(v.val, Symbol(x))
# end
#
# function create_param(x)
#     @parameters p
#     ModelingToolkit.rename(p.val, Symbol(x))
# end

to_symbol(x::Symbol) = x
to_symbol(x::AbstractString) = Symbol(x)
to_symbol(x::EzXML.Node) = Symbol(x["name"])
to_symbol(comp::Component) = comp.name

"""
    find_iv finds the unique independent variable
"""
find_iv(doc::Document) = infer_iv(doc)[:β]

function find_iv(doc::Document, comp)
    d = infer_iv(doc)
    haskey(d, to_symbol(comp)) ? d[to_symbol(comp)] : d[:β]
end

get_ivₚ(doc::Document) = create_param(find_iv(doc))

function find_components_iv(doc::Document)
    dv = Dict{Symbol, Symbol}()
    nv = Symbol[]

    for comp in components(doc)
        nodes = list_component_bvar(comp)
        ivs = unique(strip.(nodecontent.(nodes)))
        sym = to_symbol(comp)

        if length(ivs) == 1
            dv[sym] = Symbol(ivs[1])
        elseif length(ivs) == 0
            push!(nv, sym)
        else
            error("component $sym has more than two independent variables")
        end
    end
    return dv, nv
end

@memoize function infer_iv(doc::Document)
    # phase 1, find direct component ivs
    dv, dn = find_components_iv(doc)
    syms = unique(collect(values(dv)))

    if isempty(syms)
        error("no independent variable is inferred")
    end

    if isempty(dn) && length(syms) == 1
        return Dict(:β => syms[1])
    end

    groups = find_equivalence_groups(doc)
    V = ∪([groups[Var(first(x),last(x))] for x in dv]...)
    syms = unique([v.var for v in V])

    if length(syms) == 1
        return Dict(:β => syms[1])
    end

    @warn "more than one name is identified for the independent variable. It is changed to :τ."

    d = Dict(:β => :τ)
    for v in V
        d[v.comp] = v.var
    end
    return d
end

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

make_var(c, v) = Var(to_symbol(c), to_symbol(v))

names_to_varset(comp, l) = Set([make_var(comp, x) for x in l])

"""
    list_component_lhs returns a set of Var composed of all variables in component
    comp that occur on the left-hand-side of an ODE or algebraic equation
"""
function list_component_lhs(comp)
    names_to_varset(comp, find_state_names(comp)) ∪ names_to_varset(comp, find_alg_names(comp))
end

"""
    list_component_lhs returns a set of Var composed of all variables anywhere in
    the doc document that occur on the left-hand-side of an ODE or algebraic equation
"""
list_all_lhs(doc::Document) = ∪([list_component_lhs(c) for c in components(doc)]...)

"""
    find_equivalence_groups categorizes all the variables in the doc document
    based on the connections into equivalnce groups
    it returns a Dictionary of Var to groups (Set of Vars)
"""
@memoize function find_equivalence_groups(doc::Document)
    groups = Dict{Var, Set{Var}}()

    # phase 1: each variable belongs to only one set composed of the variable itself
    for c in components(doc)
        for v in list_component_variables(c)
            u = make_var(c, v)
            groups[u] = Set([u])
        end
    end

    # phase 2: equivalence groups (sets) are merged according to the connections
    for conn in connections(doc)
        for (u1,u2) in variables(conn)
            s = groups[u1] ∪ groups[u2]
            for x in s
                groups[x] = s
            end
        end
    end

    return groups
end

"""
    classify_variables divides all the variables in the doc document into
    two categorizes and returns a Dict from Var to Bool:
        true: Var is a left-hand-side variable
        false: Var is not a left-hand-side variable (should be a parameter)
"""
@memoize function classify_variables(doc::Document)
    groups = find_equivalence_groups(doc)
    lhs = list_all_lhs(doc)
    class = Dict{Var, Bool}()

    for g in groups
        n = length(last(g) ∩ lhs)
        if n > 1
            error("Equivalence group $g has more than one initial_value: $(last(g) ∩ lhs)")
        end
        class[first(g)] = (n == 1)
    end
    class
end

"""
    translate_connections translates the list of MathML connections to
    a list of ModelingToolkit equations
"""
function translate_connections(doc::Document, systems, class)
    a = []
    for conn in connections(doc)
        for (u1,u2) in variables(conn)
            if class[u1] && class[u2]
                var1 = getproperty(systems[u1.comp], u1.var)
                var2 = getproperty(systems[u2.comp], u2.var)
                push!(a, var1 ~ var2)
            end
        end
    end
    return a
end

"""
    pre_substitution generates the substitution rules to be applied to
    individual systems before structural_simplify merges them together

    it morphes variables and parameters to their correct form for ModelingToolkit
    based on the global CellML doc information
"""
function pre_substitution(doc::Document, comp, class)
    ivₚ = get_ivₚ(doc)
    ivₘ = find_iv(doc, comp)

    vars = to_symbol.(list_component_variables(comp))

    states = [create_var(x) => create_var(x, ivₚ) for x in vars if class[make_var(comp,x)]]
    params = [create_var(x) => create_param(x) for x in vars if !class[make_var(comp,x)] && !isequal(x, ivₘ)]
    ivsub =  [create_var(ivₘ) => ivₚ]

    return states ∪ params ∪ ivsub
end

function remove_rhs_diff(eqs)
    [eq.lhs => eq.rhs for eq in eqs if operation(eq.lhs) isa Differential]
end

"""
    post_substitution generates the substitution rules to be applied to
    the merged system after structural_simplify is applied

    if changes the names of the indepedent variable (iv) in each system
    to the global iv name

    TODO: this function assumes the basic iv name is the same among all systems
"""
function post_substitution(doc, systems)
    iv = find_iv(doc)
    ivₚ = get_ivₚ(doc)
    [create_param(string(s) * "₊" * string(iv)) => ivₚ for s in keys(systems)]
end

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

substitute_eqs(eqs, s) = [substitute(eq.lhs, s) ~ substitute(eq.rhs, s) for eq in eqs]

"""
    process_components is the main entry point
    it processes an doc document and returns a merged ODESystem

    use simplify=false only for debugging purposes!
"""
function process_components(doc::Document; simplify=true)
    extract_mathml.(doc.xmls)     # this is called here for the side-effect of disambiguiting equals
    infer_iv(doc)

    class = classify_variables(doc)
    systems = subsystems(doc, class)
    post_sub = post_substitution(doc, systems)

    sys = ODESystem(
        translate_connections(doc, systems, class),
        get_ivₚ(doc),
        systems=collect(values(systems))
    )

    if simplify
        sys = structural_simplify(sys)
        @set! sys.eqs = substitute_eqs(equations(sys), post_sub)
    end

    return sys
end

"""
    subsystems returns a Dict of the subsystems (keys are Symbols), each one
    based on one CellML component of the doc document.

    class is the output of classify_variables
"""
function subsystems(doc::Document, class)
    Dict{Symbol,ODESystem}(
        to_symbol(comp) => process_component(doc, comp, class)
        for comp in components(doc)
    )
end

"""
    process_component converts a single CellML component to an ODESystem

    comp in the name of the component
    class is the output of classify_variables
"""
function process_component(doc::Document, comp, class)
    math = list_component_math(comp)
    pre_sub = pre_substitution(doc, comp, class)

    if length(math) == 0
        eqs = Equation[]
    else
        eqs = vcat(parse_node.(math)...)
        eqs = substitute_eqs(eqs, pre_sub)
        sub_rhs = remove_rhs_diff(eqs)
        eqs = [eq.lhs ~ substitute(eq.rhs, sub_rhs) for eq in eqs]
    end

    ivₚ = get_ivₚ(doc)
    ps = [last(x) for x in values(pre_sub) if last(x) isa Sym && !isequal(last(x), ivₚ)]
    states = [last(x) for x in values(pre_sub) if !(last(x) isa Sym)]

    ODESystem(eqs, ivₚ, states, ps; name=to_symbol(comp))
end

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

function collect_initiated_values(doc::Document)
    vars = Dict{Var, Float64}()
    syms = Dict{Var, Var}()

    for comp in components(doc)
        for v in list_initiated_variables(comp)
            val = Meta.parse(v["initial_value"])
            if val isa Symbol
                syms[make_var(to_symbol(comp), v)] = make_var(to_symbol(comp), val)
            else
                vars[make_var(to_symbol(comp), v)] = Float64(val)
            end
        end
    end
    vars, syms
end

function split_sym(sym)
    # # remove dependency, e.g., (time) from the end
    s = replace(string(sym), r"\([^\)]*\)" => "")
    make_var(split(s, "₊")...)
end

find_sys_p(doc::Document, sys) = find_list_value(doc, parameters(sys))
find_sys_u0(doc::Document, sys) = find_list_value(doc, states(sys))

function find_list_value(doc::Document, names)
    vars, syms = collect_initiated_values(doc)
    varkeys = Set(keys(vars))
    groups = find_equivalence_groups(doc)

    vals = []

    for x in names
        u = split_sym(x)
        u = haskey(syms,u) ? syms[u] : u
        var = groups[u] ∩ varkeys
        if length(var) == 0
            error("value of $x is not found")
        end
        push!(vals, x => vars[first(var)])

    end
    return vals
end


text:
using VLMinorityGameKit
using BSON

# setup call -
memory = 10
numberOfStrategiesPerAgent = 5
numberOfAgents = 1001

# setup -
path_to_game_world = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/GW-M10-S5-NA1001.bson"
# path_to_game_world = "/Users/jeffreyvarner/Desktop/proposals/VLMinorityGameKit.jl/test/data/GW-M10-S5-NA5.bson"

# build the game world -
# args: numberOfAgents::Int64, agentMemorySize::Int64, numberOfStrategiesPerAgent::Int64; actions::Array{Int64,1} = [-1,0,1]
game_world = build_minority_game_world(VLBasicMinorityGameWorld; numberOfAgents, memory, numberOfStrategiesPerAgent)

# write the gameworld to disk -
bson(path_to_game_world, gameworld=game_world)