text:
module Day23

export get_inputs, get_solution1, get_solution2

using CircularArrays: CircularArray
using ProgressMeter: @showprogress

## Input getting
function get_inputs()
    val = read_input(389125467)
    test_input1 = [(val, 10), (val, 100)]
    test_output1 = [92658374, 67384529]
    test_input2 = val
    test_output2 = nothing
    data = (read_input(joinpath(@__DIR__, "input.txt")), 100)
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

read_input(n::Number) = n |> digits |> reverse!
read_input(io) = parse(Int, only(readlines(io))) |> read_input

## Solution functions
function permutation_indices!(perm_vector, cups)
    for (i, cup) in enumerate(cups)
        perm_vector[cup] = i
    end
    return perm_vector
end

function insertall!(a, i, vals)
    for val in reverse!(vals)
        insert!(a, i, val)
    end
    return a
end

function play_game(cups, nturns)
    cups = CircularArray(copy(cups))
    ncups = length(cups)
    indices = sortperm(cups.data)
    last_cup = last(cups)
    pickup_cups = zeros(Int, 3)
    last_index = 0

    @showprogress 1 for turn in 1:nturns
        # last_index = findfirst(==(last_cup), cups)

        current_cup = cups[mod(last_index+1, ncups)]

        # println("-- move $turn --")
        # println("cups:    $cups")
        # println("indices: $indices")
        # println("current cup: $current_cup")

        pickup_indices = last_index+2:last_index+4
        pickup_index = first(pickup_indices)
        pickup_cups .= @view cups[pickup_indices]

        destination_cup = current_cup
        for i in 1:4
            d = mod1(destination_cup-i, ncups)
            if d ∉ pickup_cups
                destination_cup = d
                break
            end
        end

        # println("pick up: $pickup_cups")
        # println("destination: $destination_cup\n")

        destination_index = findfirst(==(destination_cup), cups)
        # destination_index = indices[destination_cup]
        if destination_index < pickup_index
            destination_index += ncups
        end

        replace_indices = pickup_index:destination_index-1
        @inbounds for index in replace_indices
            if index > ncups
                index = mod1(index, ncups)
            end
            cup = cups[index+3]
            cups.data[index] = cup

            indices[cup] = mod1(index, ncups)
        end
        # cups[pickup_index:destination_index-1] .= @view cups[pickup_index+3:destination_index+2]

        putdown_indices = mod1.(destination_index-2:destination_index, ncups)
        cups[putdown_indices] .= pickup_cups
        indices[pickup_cups] .= putdown_indices

        last_cup = current_cup
        last_index = findfirst(==(last_cup), cups)
        # last_index = indices[last_cup]
        # last_index += 1
        # if (last_index+1) in (pickup_index:destination_index+3)
        #     last_index -= 3
        # end
    end

    return cups
end

# function play_game(cups, turns, ncups=9)
#     cups = copy(cups)
#     perm_vector = similar(cups)
#     permutation_indices!(perm_vector, cups)
#     last_cup = last(cups)

#     for turn in 1:turns
#         last_i = findfirst(==(last_cup), cups)+1
#         # println("-- move $turn --")
#         # println("cups: $cups")
#         i = mod1.(last_i:last_i+3, ncups)
#         current_cup = cups[i[1]]
#         last_cup = current_cup
#         three_inds = i[2:4]
#         three_cups = cups[three_inds]
#         foreach(idx-> deleteat!(cups, idx), sort!(three_inds, rev=true))
#         destination_cup = current_cup
#         for i in 1:4
#             destination_cup = mod1(destination_cup-1, ncups)
#             destination_cup ∉ three_cups && break
#         end
#         # println("pick up: $three_cups")
#         # println("destination: $destination_cup\n")
#         permutation_indices!(perm_vector, cups)
#         i_dest = perm_vector[destination_cup]
#         # i_dest = findfirst(==(destination_cup), cups)
#         insertall!(cups, mod1(i_dest+1, ncups), three_cups)
#         # permutation_indices!(perm_vector, cups)
#     end

#     return cups
# end

function get_solution(cups, nturns)
    cups = play_game(cups, nturns)
    i_one = findfirst(==(1), cups)
    inds = mod1.(i_one+1:i_one+8, length(cups))
    return parse(Int, join(cups[inds]))
end

# Part 1
get_solution1(data::Vector{<:Tuple}) = get_solution1.(data)
get_solution1(data::Tuple) = get_solution1(data...)
get_solution1(cups, turns) = get_solution(cups, turns)

# Part 2
get_solution2(data::Tuple) = get_solution2(data[1])
function get_solution2(cups)
    cups = [cups; (length(cups)+1):1_000_000]
    cups = play_game(cups, 10_000_000)
    i = findfirst(==(1), cups)
    return prod(cups[i+1:i+2])
end

end


text:
# only linear PHS for now
# coupling using gyrator-type interconnection
function coupled_gyrator(ph1 :: Phs, ports1, ph2 :: Phs, ports2, couple_matrix)
	N1 = size(ph1.J,1); # number of states
	N2 = size(ph2.J,1);

	# split the input ports that are not used for interconnection
	# and creates the new B matrix:
	inputs1 = [];
	for input = 1:size(ph1.B,2)
		if ~(input in ports1)
			inputs1 = push!(inputs1,input);
		end
	end
	inputs2 = [];
	for input = 1:size(ph2.B,2)
		if ~(input in ports2)
			inputs2 = push!(inputs2,input);
		end
	end
	Bnew = blkdiag(ph1.B[:,inputs1],ph2.B[:,inputs2]);
	Bnew[1:N1,(length(inputs1)+1):end] = ph1.B[:,ports1]*couple_matrix*ph2.D[ports2,inputs2];
	Bnew[(N1+1):end, 1:length(inputs1)] = -ph2.B[:,ports2]*couple_matrix*ph1.D[ports1,inputs1] ;

	Jnew = blkdiag(ph1.J, ph2.J);
	Jnew[1:N1, (N1+1):end] = ph1.B[:,ports1]*couple_matrix*(transpose(ph2.B[:,ports2]))
	Jnew[(N1+1):end, 1:(N1)] = -ph2.B[:,ports2]*couple_matrix*(transpose(ph1.B[:,ports1]))

	Dnew = zeros(length(inputs1)+length(inputs2),length(inputs1)+length(inputs2));
	Dnew[1:length(inputs1), 1:length(inputs1)] = ph1.D[inputs1,inputs1]
	Dnew[length(inputs1)+1:end, length(inputs1)+1:end] = ph2.D[inputs2,inputs2]
	Dnew[1:length(inputs1), length(inputs1)+1:end] = ph1.D[inputs1,ports1]*couple_matrix*ph2.D[ports2, inputs2];
	Dnew[length(inputs1)+1:end, 1:length(inputs1)] = - ph2.D[inputs2, ports2]*couple_matrix*ph1.D[ports1,inputs1];
	Qnew = blkdiag((ph1.Q), (ph2.Q));
	phsnew = Phs(Jnew, Bnew, Dnew, Qnew)
	phsnew.TransfMatrix = blkdiag(ph1.TransfMatrix,ph2.TransfMatrix)

	# preserves the names of the states:
	N1 = size(ph1.TransfMatrix,1)
	for key in keys(ph2.StatesNames)
		ph2.StatesNames[key] += N1
	end
	phsnew.StatesNames = merge(ph1.StatesNames, ph2.StatesNames)

	# remove the names of the coupled inputs
	inp1 = copy(ph1.InputsNames)
	deleteat!(inp1, ports1);
	inp2 = copy(ph2.InputsNames)
	deleteat!(inp2, ports2);
	phsnew.InputsNames = [inp1; inp2]

	phsnew.R = blkdiag(ph1.R, ph2.R)
	return phsnew
end

function coupled_transformer(ph1 :: Phs, ports1, ph2 :: Phs, ports2, couple_matrix)
	N1 = size(ph1.J,1); # number of states
	N2 = size(ph2.J,1);

	# split the input ports that are not used for interconnection
	# and creates the new B matrix:
	inputs1 = [];
	for input = 1:size(ph1.B,2)
		if ~(input in ports1)
			inputs1 = push!(inputs1,input);
		end
	end
	inputs2 = [];
	for input = 1:size(ph2.B,2)
		if ~(input in ports2)
			inputs2 = push!(inputs2,input);
		end
	end
	Bnew = blkdiag(ph1.B[:,inputs1],ph2.B[:,inputs2]);
	#Bnew[1:N1,(length(inputs1)+1):end] = ph1.B[:,ports1]*couple_matrix*ph2.D[ports2,inputs2];
	#Bnew[(N1+1):end, 1:length(inputs1)] = -ph2.B[:,ports2]*couple_matrix*ph1.D[ports1,inputs1] ;

	Jnew = blkdiag(ph1.J, ph2.J);
	#Jnew[1:N1, (N1+1):end] = ph1.B[:,ports1]*couple_matrix*(ph2.B[:,ports2].')
	#Jnew[(N1+1):end, 1:(N1)] = -ph2.B[:,ports2]*couple_matrix*(ph1.B[:,ports1].')

	#Dnew = zeros(length(inputs1)+length(inputs2),length(inputs1)+length(inputs2));
	#Dnew[1:length(inputs1), length(inputs1)+1:end] = ph1.D[inputs1,ports1]*couple_matrix*ph2.D[ports2, inputs2];
	#Dnew[length(inputs1)+1:end, 1:length(inputs1)] = - ph2.D[inputs2, ports2]*couple_matrix*ph1.D[ports1,inputs1];
	G = [ph1.B[:,ports1] ; ph2.B[:,ports2]*couple_matrix]
	GD = ph1.D[ports1,ports1] - ph2.D[ports2,ports2]*couple_matrix
	Dnew = blkdiag(ph1.D[inputs1,inputs1], ph2.D[inputs2,inputs2])
	Qnew = blkdiag((ph1.Q), (ph2.Q));
	phnew = Phs(Jnew, Bnew, Dnew, Qnew)
	set_constraint!(phnew, G, GD)
	phnew.TransfMatrix = blkdiag(ph1.TransfMatrix,ph2.TransfMatrix)

	# preserves the names of the states:
	N1 = size(ph1.TransfMatrix,1)
	for key in keys(ph2.StatesNames)
		ph2.StatesNames[key] += N1
	end
	phnew.StatesNames = merge(ph1.StatesNames, ph2.StatesNames)
	phnew.R = blkdiag(ph1.R, ph2.R)
	return phnew
end

ext: jl

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
using BlockCartesianGrids
using Test

max_stars_repo_name: lcw/BlockCartesianGrids

max_issues_repo_name: lcw/BlockCartesianGrids

max_forks_repo_name: lcw/BlockCartesianGrids


text:
# %% Sensititivy to trajectory length
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

N = Int(1e5)
Ds = range(4; step = 2, length = 6)

qH = 1
qC = 2
data = :lorenz96_chaotic

# Calculate values for H, C
eHs, eCs, Hs, Cs = [Vector{Float64}[] for i in 1:4]
for i in 1:length(Ds)
    # These are the parameters that change between the different elements of the plot
    D = Ds[i]

    # Here we simply pack all parameters into a dictionary
    # (other parameters are (probably) globals)
    params = @strdict N D qH qC data

    # This is the main call that calculates everything
    output, s = produce_or_load(
        datadir("main"), params, make_C_H;
        prefix = string(data), suffix = "jld2", force = true,
        ignores = ["data"],
    )
    @unpack eH, eC, H, C = output
    push!(eHs, eH); push!(Hs, H); push!(eCs, eC); push!(Cs, C)
end

# Do the actual plot
labels = ["\$D=$(D)\$" for D in Ds]
legendtitle = "impact of dimensionality \$D\$ (Lorenz96 system)"

fig, axs = mainplot(
    Hs, Cs, es, labels, legendtitle;
    qH, qC, tol = 0.25,
    # dimension_fit = FractalDimension.logarithmic_corrected_fit_lsqfit
)

# wsave(plotsdir("paper", "dimensionality"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
# sample sequences of ergms with different parameters' values from R package ergm
# and test the PseudoLikelihoodScoreDrivenERGM filter
using Utilities,AReg,StaticNets,DynNets , JLD,MLBase,StatsBase,CSV, RCall
 using PyCall; pygui(); using PyPlot
 using ErgmRcall

# load packages needed in R
R"rm(list = ls())
   library(statnet)
    library(ergm)
    library(sna)
    library(coda)
    library(network)
    sessionInfo()"

tot_Nsample = 200
N_each_chunk = 10
#divide the simulations in chuncks
for n_chunk = 1:Int(tot_Nsample/N_each_chunk)
T = 110
N = 100
UM = [-3 , 2.5]
B_Re  = 0.95
A_Re  = 0.05
indTvPar = trues(2)
indTargPar = trues(2)
Nsample = N_each_chunk
estPar_all = fill(fill(ones(2),2),Nsample)
obsT_all = fill(fill(ones(2,2),T),Nsample)

convFlag  = falses(Nsample)
n=1
while n <= Nsample
    # sample the dgp
    print(n)

    R"rm(list = ls())
       library(statnet)
        library(ergm)
        library(sna)
        library(coda)
        library(network)
        sessionInfo()"
    try
        obsT,statsT,fVecT = ErgmRcall.sdergmDGPpaper_Rcall(  DynNets.fooSdErgmDirBinGlobalPseudo, indTvPar ,T,N, UM = UM, B_Re = B_Re,A_Re = A_Re)
        model = DynNets.SdErgmDirBinGlobalPseudo(obsT,DynNets.fooGasPar,indTvPar,"")
        estPar,convFlag[n],UM_est,startPoint = DynNets.estimate(model;indTvPar = indTvPar, indTargPar = indTargPar)
        obsT_all[n] = obsT
        estPar_all[n] = estPar
        if !isnan(estPar[2][1])
            n = n+1
        end
    catch
        # clear all R variables
    end
end
# Save all Sampled Data
save_fold = "./data/estimatesTest/sdergmTest/distrib_static_pars/"
 @save(save_fold*"sdergm_dgp__$(N)_T_$(T)_Sample_$(Nsample)_UM_$(UM)_intTV_$(Int.(indTvPar[:]))_intTarg_$(Int.(indTargPar[:]))_N_each_chunk_$(N_each_chunk)_chunk_$(n_chunk).jld",
              obsT_all, estPar_all, UM,B_Re,A_Re,convFlag, Nsample,T,N, indTvPar, indTargPar )

end
#
#
# using Utilities,AReg,StaticNets,DynNets , JLD,MLBase,StatsBase,CSV, RCall
#  using PyCall; pygui(); using PyPlot
#  using JLD2,Utilities, GLM
#  ## load R MCMC simulation and estimates and estimate sdergmTest
#
#  if false
#  Nsample = 100
#  dgpType = "sin"
#  T = 50
#  N = 50
#  Nterms = 2
#  Nsteps1 ,Nsteps2 = 0,1
#  load_fold = "./data/estimatesTest/sdergmTest/R_MCMC_estimates/"
#  @load(load_fold*"test_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE.jld",
#              stats_T, changeStats_T,estParSS_T,sampledMat_T ,parDgpT,Nsample)
#  end
#
#  onlyTest = false
#  targetAllTv = true
#  tvParFromTest = false; pValTh = 0.05
#  gasParEst = fill(fill(Float64[],2), Nsample)
#  startPointEst = fill(fill(0.0,2), Nsample)
#  staticEst = fill(fill(0.0,2), Nsample)
#  filtPar_T_Nsample = zeros(Nterms,T,Nsample)
#  pVals_Nsample = zeros(Nterms,Nsample)
#  scoreAutoc_Nsample = zeros(Nterms,Nsample)
#  convFlag = falses(Nsample)
#  tmp = Array{Array{Float64,2},2}(T,Nsample); for t=1:T,s=1:Nsample tmp[t,s] =  changeStats_T[t][s];end;#changeStats_T = tmp
#  for n=1:Nsample
#      @show(n)
#       pVals_Nsample[:,n],tmpInfo,staticEst[n] =
#          DynNets.pValStatic_SDERGM( DynNets.SdErgmDirBinGlobalPseudo(tmp[:,n],DynNets.fooGasPar,falses(Nterms),""))
#       @show(tmpInfo)
#      if tvParFromTest
#          indTvPar =  pVals_Nsample[:,n] .< pValTh
#      else
#          indTvPar = BitArray([true,true])
#      end
#      if targetAllTv
#          indTargPar = indTvPar
#      else
#           indTargPar =BitArray([false,false])
#      end
#     model = DynNets.SdErgmDirBinGlobalPseudo(tmp[:,n],DynNets.fooGasPar,indTvPar,"")
#
#     if .!onlyTest
#     estPar,convFlag[n],UM,startPoint = DynNets.estimate(model;UM =meanSq(estParSS_T[n,:,:],1),indTvPar = indTvPar, indTargPar = indTargPar)
#
#
#     gasParEst[n] = estPar # store gas parameters
#     startPointEst[n] = startPoint
#     gasParVec = zeros(sum(indTvPar)*3); for i=0:(sum(indTvPar)-1) gasParVec[1+3i : 3(i+1)] = estPar[indTvPar][i+1]; end
#     constParVec = zeros(sum(!indTvPar)); for i=1:sum(.!indTvPar) constParVec[i] = estPar[.!indTvPar][i][1]; end
#     gasFiltPar , pseudolike = DynNets.score_driven_filter_or_dgp(model,gasParVec,indTvPar;
#                               vConstPar = constParVec,ftot_0= startPoint )
#
#
#     filtPar_T_Nsample[:,:,n] = gasFiltPar
#     end
#  end
#
#  save_fold = "./data/estimatesTest/sdergmTest/"*
#             "gas_MCMC_comparison_estimates/"
#  @save(save_fold*"test_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE_target_$(targetAllTv).jld" ,
#      stats_T, changeStats_T,estParSS_T,sampledMat_T ,parDgpT,
#      Nsample,T,N,filtPar_T_Nsample,gasParEst,convFlag,pVals_Nsample,scoreAutoc_Nsample,staticEst)
#
#
# if false
#
#     sum(convFlag)
#     figure()
#      legTex = ["Constant ERGM";"SDERGM"; "Sequence of ERGM"; "DGP"]
#      for n=1:Nsample
#      indTvPar = BitArray([true,true])
#      model = SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,indTvPar,"")
#      estPar = gasParEst[n] # store gas parameters
#      gasFiltPar  = filtPar_T_Nsample[:,:,n]
#      parInd = 1
#      subplot(1,2,1);plot(1:T,ones(T)*staticEst[n][parInd],"-b")
#                     plot(1:T,gasFiltPar[parInd,:],"r")
#                     plot(1:T,estParSS_T[n,:,parInd],".b")
#                     plot(1:T,parDgpT[parInd,:],"k",linewidth=5)
#
#      parInd = 2
#      subplot(1,2,2);plot(1:T,ones(T)*staticEst[n][parInd],"-b")
#                     plot(1:T,gasFiltPar[parInd,:],"r")
#                     plot(1:T,estParSS_T[n,:,parInd],".b")
#                     plot(1:T,parDgpT[parInd,:],"k",linewidth=5)
#      end
#      namePar1 = "Number of Links"
#      namePar2 = "GWESP"
#      subplot(1,2,1);    title(namePar1); legend(legTex)
#      subplot(1,2,2);  title(namePar2 ); legend(legTex)
#
#     #Plotta info relative ai tests
#     figure()
#          subplot(2,2,1); parInd = 1;   plt[:hist]((pVals_Nsample[parInd,:]), bins=logspace(minimum(log10(pVals_Nsample[1,:])),maximum(log10(pVals_Nsample[parInd,:])), 20))
#          if maximum((pVals_Nsample[parInd,:])) >pValTh
#               axvspan((pValTh),maximum((pVals_Nsample[parInd,:])),color = "r",alpha = 0.1);
#          end
#               xscale("log");
#          subplot(2,2,2);  parInd = 2;   plt[:hist]((pVals_Nsample[parInd,:]), bins=logspace(minimum(log10(pVals_Nsample[1,:])),maximum(log10(pVals_Nsample[parInd,:])), 20))
#          if maximum((pVals_Nsample[parInd,:])) >pValTh
#               axvspan((pValTh),maximum((pVals_Nsample[parInd,:])),color = "r",alpha = 0.1);
#          end
#               xscale("log");
#          subplot(2,2,3); plot((pVals_Nsample[1,:]),(pVals_Nsample[2,:]),".")
#          axvline((pValTh));axhline((pValTh))
#          xscale("log");yscale("log")
#          xlabel("P-Value " * namePar1)
#          ylabel("P-Value " * namePar2)
# end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# Vector of (eclassid, position_of_literal_in_eclass_nodes)
const Sub = Tuple{Vector{Tuple{Int64, Int64}}, Dict{Any, Type}}

haseclassid(sub::Sub, p::PatVar) = first(sub[1][p.idx]) >= 0
geteclassid(sub::Sub, p::PatVar) = first(sub[1][p.idx])

hasliteral(sub::Sub, p::PatVar) = last(sub[1][p.idx]) > 0
getliteral(sub::Sub, p::PatVar) = last(sub[1][p.idx])

hastermtype(sub::Sub, p::Any) = haskey(sub[2], p)
gettermtype(sub::Sub, p::Any) = sub[2][p]

## ====================== Instantiation =======================

function instantiate(g::EGraph, pat::PatVar, sub::Sub, rule::Rule)
    if haseclassid(sub, pat)
        ec = geteclass(g, geteclassid(sub, pat))
        if hasliteral(sub, pat)
            pos = getliteral(sub, pat)
            node = ec.nodes[pos]
            @assert arity(node) == 0
            return node.head
        end
        return ec
    else
        error("unbound pattern variable $pat in rule $rule")
    end
end

function instantiate(g::EGraph, pat::PatLiteral{T}, sub::Sub, rule::Rule) where T
    pat.val
end

function instantiate(g::EGraph, pat::PatTypeAssertion, sub::Sub, rule::Rule)
    instantiate(g, pat.name, sub, rule)
end

# # TODO CUSTOMTYPES document how to for custom types
function instantiateterm(g::EGraph, pat::PatTerm,  T::Type{Expr}, sub::Sub, rule::Rule)
    Expr(pat.head, map(x -> instantiate(g, x, sub, rule), pat.args)...)
end

# FIXME
function instantiate(g::EGraph, pat::PatTerm, sub::Sub, rule::Rule)
    # println(sub)
    # for (pp, tt) ∈ sub.termtypes
    #     println("$pp $(hash(pp))")
    #     println("$pat $(hash(pat))")
    #     println(pp == pat)
    # end
    if hastermtype(sub, pat.head)
        T = gettermtype(sub, pat.head)
        # println(T)
        instantiateterm(g, pat, T, sub, rule)
    else
        # dump(pat)
        # println("$pat has no type in $sub")
        instantiateterm(g, pat, Expr, sub, rule)
    end
end


text:
using Test
using MyTestPkg

my_tests = [
    "test_sum"
]

@testset "All tests" begin
    for test in my_tests
        println(test)
        include("$(test).jl")
    end
end


text:
converter = AvgCnvFixedDC(690.0, 2750000.0) #S_rated goes in Watts
dc_source = FixedDCSource(600.0) #Not in the original data, guessed.
filter = LCLFilter(0.08, 0.003, 0.074, 0.2, 0.01)
pll = PLL(500.0, 0.084, 4.69)

virtual_H = VirtualInertia(2.0, 400.0, 20.0, 2*pi*50.0)
Q_control = ReactivePowerDroop(0.2, 1000.0)
outer_control = VirtualInertiaQdroop(virtual_H, Q_control)

vsc = CurrentControl(1.27, 14.3, 0.59, 736.0, 50.0, 0.5, 0.2,
0.00)

DynInverter(1,
            :DARCO,
            nodes_OMIB[2],
            1.0,
            1.0,
            0.4,
            0.0,
            2.75e3,
            converter,
            outer_control,
            vsc,
            dc_source,
            pll,
            filter)


text:
### Julia module for the Compressed Stack structure ###
__precompile__() # Precompile option for the module. Is recompiled after changes
module CompressedStacks

## Import/Export
export CompressedStack, NormalStack, run!,runCountingPop!, top, top_explicit

## Basic types and constructors for CompressedStack
include("base.jl")
## NormalStack
include("normal.jl")
## I/O for Stacks structure
include("io.jl")

## Action for Compressed Stack's internals mechanisms
# Basic access, test emptiness, update blocks
include("access.jl")
# Pushes and pops
include("push.jl")
include("pop.jl") # reconstruct.jl is nested in
# Run the stacks
include("run.jl")

## Specialized types of Compressed Stacks

end ### of CompressedStacks module


text:
# rand
import Base.rand

typealias RandTypes{T,C} Union{AbstractRGB{T},
                               TransparentRGB{C,T},
                               AbstractGray{T},
                               TransparentGray{C,T}}
typealias RandTypesFloat{C,T<:AbstractFloat} RandTypes{T,C}

rand{G<:AbstractGray}(::Type{G}) = G(rand())
rand{G<:TransparentGray}(::Type{G}) = G(rand(), rand())
rand{C<:AbstractRGB}(::Type{C}) = C(rand(), rand(), rand())
rand{C<:TransparentRGB}(::Type{C}) = C(rand(), rand(), rand(), rand())
function rand{C<:RandTypesFloat}(::Type{C}, sz::Dims)
    reinterpret(C, rand(eltype(C), (sizeof(C)÷sizeof(eltype(C)), sz...)), sz)
end
function rand{C<:RandTypes{U8}}(::Type{C}, sz::Dims)
    reinterpret(C, rand(UInt8, (sizeof(C), sz...)), sz)
end
function rand{C<:RandTypes{U16}}(::Type{C}, sz::Dims)
    reinterpret(C, rand(UInt16, (sizeof(C)>>1, sz...)), sz)
end

# Mapping a function over color channels
"""
    mapc(f, rgb) -> rgbf
    mapc(f, rgb1, rgb2) -> rgbf

`mapc` applies the function `f` to each color channel of the input
color(s), returning an output color in the same colorspace.

# Examples:

    julia> mapc(x->clamp(x,0,1), RGB(-0.2,0.3,1.2))
    RGB{Float64}(0.0,0.3,1.0)

    julia> mapc(max, RGB(0.1,0.8,0.3), RGB(0.5,0.5,0.5))
    RGB{Float64}(0.5,0.8,0.5)

    julia> mapc(+, RGB(0.1,0.8,0.3), RGB(0.5,0.5,0.5))
    RGB{Float64}(0.6,1.3,0.8)
"""
mapc{C<:AbstractGray}(f, c::C) = base_color_type(C)(f(gray(c)))
mapc{C<:TransparentGray}(f, c::C) = base_colorant_type(C)(f(gray(c)), f(alpha(c)))
mapc{C<:Color3}(f, c::C) = base_color_type(C)(f(comp1(c)), f(comp2(c)), f(comp3(c)))
mapc{C<:Transparent3}(f, c::C) = base_colorant_type(C)(f(comp1(c)), f(comp2(c)), f(comp3(c)), f(alpha(c)))

mapc(f, x, y) = _mapc(_same_colorspace(x,y), f, x, y)
_mapc{C<:AbstractGray}(::Type{C}, f, x, y) = C(f(gray(x), gray(y)))
_mapc{C<:TransparentGray}(::Type{C}, f, x, y) = C(f(gray(x), gray(y)), f(alpha(x), alpha(y)))
_mapc{C<:Color3}(::Type{C}, f, x, y) = C(f(comp1(x), comp1(y)), f(comp2(x), comp2(y)), f(comp3(x), comp3(y)))
_mapc{C<:TransparentRGB }(::Type{C}, f, x, y) = C(f(red(x), red(y)), f(comp2(x), comp2(y)), f(comp3(x), comp3(y)), f(alpha(x), alpha(y)))

_same_colorspace(x::Colorant, y::Colorant) = _same_colorspace(base_colorant_type(x),
                                                              base_colorant_type(y))
_same_colorspace{C<:Colorant}(::Type{C}, ::Type{C}) = C
@noinline _same_colorspace{C1<:Colorant,C2<:Colorant}(::Type{C1}, ::Type{C2}) =
    throw(ArgumentError("$C1 and $C2 are from different colorspaces"))


text:
#
"""
 FWD Pass: NonLin layer -> Lin Layer
 BWD Pass: GD on NonLin layer, Lin layer

"""

using Revise
using diffMesh
import Zygote
using NNlib, Flux, LinearAlgebra, Plots
#--------------------------------------#
np  = 2001 # data points
ng  = 31   # FEM nodes

function udata(x)
    u   = 2 .+ zero(x)
    u .+= @. sin(1*pi*x) * (-1 <= x <= 0)
    u .+= @. sin(4*pi*x) * ( 0 <  x <= 1)
    return u
end
#--------------------------------------#
function addEndPts(xe::AbstractArray
                  ,x0::Number,x1::Number)

    n = length(xe) + 2
    A,R = Zygote.ignore() do
        A = zeros(n,2)
        A[1,1] = 1
        A[n,2] = 1
        R = restrict(n)
        return A,R
    end

    xg = A * [x0,x1] + R'*xe

    return xg
end
#--------------------------------------#
function model(xe,ug,x)

    xg = addEndPts(xe,-1.0,1.0)

    return reluHat(xg,ug,x)
end
#--------------------------------------#
function llsq(xg)

    NN = reluHatFunc(xg)
    A  = NN(xp')'
    linOp(u) = reluHat(NN,u,xp)

#   ug = udata(xg) # wrong   least sq. for mse loss minimization
    ug = A \ up    # correct least sq. for mse loss minimization

    return ug
end
#--------------------------------------#
function loss()
    fp   = model(xe,ug,xp)
    loss = Flux.Losses.mse(fp,up)
    return loss
end
#--------------------------------------#
function makeplot()
    fp = model(xe,ug,xp)

    xg = addEndPts(xe,-1,1)
    NN = reluHatFunc(xg)
#   ug;= llsq(xg) # defined in global scope

    plt = plot()
    plt = plot!(title="Gradient Descent",xlabel="x",ylabel="y")
    plt = plot!(ylims=(-1,4))

    plt = plot!(xp,up,width=3,label="data")
    plt = plot!(xp,fp,width=3,label="prediction")
#   plt = scatter!(xg,ug,width=4,color="red",label=false)

    # hat functions
    NNp = NN(xp')
    for i=1:ng
        scale = ug[i]
#       scale = 1.0
        global plt = plot!(xp,scale .* NNp[i,:]
                          ,width=2,color="black",label=false)
    end

    #savefig(plt,"diffMesh.png")
    display(plt)
    return plt
end
#--------------------------------------#
function callback(;doplot=true)

    l2 = loss()
    fp = model(xe,ug,xp)
    er = norm(up .- fp,Inf)
    println("Error ∞ norm: $er, loss: $l2")

    if(doplot) makeplot() end
    return false
end
#--------------------------------------#

# data
xp = Array(range(-1,stop=1,length=np))
up = udata(xp)

# initialize
xg = sort(2*rand(ng).-1)
ug = rand(length(xg))
#xg = Array(range(-1,stop=1,length=ng))
#ug = udata(xg)

xe = xg[2:end-1] # w/o end points
fp = model(xe,ug,up)

callback(doplot=true)

ps = Flux.params(xe,ug)
data = Iterators.repeated((),100)
opt = ADAM(1e-2)

gs = gradient(()->loss(),ps)

#anim = Animation()
for d in data
    global gs .= gradient(()->loss(),ps)
    Flux.Optimise.update!(opt,ps,gs)
    callback(doplot=true)
    #frame(anim)
end
#gif(anim,"examples/GD.gif",fps=30)
#--------------------------------------#
nothing


text:
queue = PriorityQueue{String,Priority}()

enqueue!(queue, "1"=>Priority(2.0))
enqueue!(queue, "2"=>Priority(2.0))
enqueue!(queue, "3"=>Priority(4.0))

x = dequeue!(queue)

@test (x== "1" || x == "2")
y = dequeue!(queue)

@test (y== "1" || y == "2")
@test dequeue!(queue) == "3"

enqueue!(queue, "4", Priority(5,1))
enqueue!(queue, "5", Priority(5,0))

@test dequeue!(queue) == "5"
@test dequeue!(queue) == "4"


text:
module MPFR_wrap

using Base.MPFR: MPFR, MPFRRoundingMode
using Base.MPFR: DEFAULT_PRECISION, ROUNDING_MODE
using Base.MPFR: CulongMax, ClongMax, CdoubleMax

include("assignment.jl")
include("basic_arithmetic.jl")
include("comparison.jl")
include("special_functions.jl")
include("integer_reminder.jl")

end # module


text:
const mac_gdrive = "/Volumes/GoogleDrive/My Drive"
const linux_gdrive = "/home/beatriz/mainen.flipping.5ht@gmail.com"
const files_dir = "Flipping/Datasets/FMRIOpto"
const figs_dir = "Flipping/Datasets/FMRIOpto/figures/"
## load
if ispath(linux_gdrive)
    ongoing_dir = linux_gdrive
else
    ongoing_dir = mac_gdrive
end

files_loc = joinpath(ongoing_dir,files_dir)
figs_loc = joinpath(ongoing_dir,figs_dir)

const columns_types = Dict(
    :Mouse => String,
    :Genotype => String,
    :Area => String,
    :Concentration => Float64,
    :Stim => Bool,
    :Repetition_nr => Int,
    :Trial_nr => Int,
    :Block => Int,
    :BOLD_90 => Float64,
    :BOLD_85 => Float64,
    :BOLD_0 => Float64,
    :Frame => Int,
    :Valid_frs => Bool,
    :Outlier_count_per_rep => Int
    )
using FMRIOpto
fulldata =  CSV.read(joinpath(files_loc,"odour_dose_data.csv"), DataFrame;
    types = columns_types,
    missingstring = "NaN")
mapcols!(Array, fulldata)
fulldata.Time = fulldata.Frame .+ 3
## average response per concentraition and stim over all presentation
gd0 = groupby(fulldata,[:Genotype, :Mouse, :Stim,:Area, :Concentration, :Time])
df0 = combine(gd0, [:BOLD_90, :BOLD_85,:BOLD_0] .=> mean .=> [:BOLD_90, :BOLD_85,:BOLD_0])

## define bold
fdata0 = filter(r -> r.Genotype == "sert" &&
    !(r.Area in ["pPIR", "OT", "NDB", "LHA", "COA"]) &&
    r.Concentration != 0.0001,
    # r.Mouse != "FM14",
    df0)
gd1 = groupby(fdata0, [:Mouse, :Area, :Stim, :Concentration])
df1 = combine(gd1, :BOLD_90 => mean => :BOLD_90)
plt, res, m, mOFC, mOB = analyze_bold(df1,:BOLD_90)
plt
res
m
mOFC
mOB
##
df2 = combine(groupby(fdata0, [:Genotype, :Mouse, :Stim,:Area, :Concentration]),
    [:BOLD_90, :BOLD_85,:BOLD_0] .=> max_bold .=> [:BOLD_90, :BOLD_85,:BOLD_0])
plt, res, m, mOFC, mOB, maPIR = analyze_bold(df2,:BOLD_90)
plt
res
m
mOFC
mOB


text:
import Base.sum
function sum(A::Assoc, dim::Int64)
    if(dim == 1)
        return condense(Assoc(promote([1.0],A.row)[1],A.col, promote([1.0],A.val)[1],sparse(sum(A.A,1))))
    elseif (dim == 2)

        return condense(Assoc(A.row,promote([1.0],A.col)[1], promote([1.0],A.val)[1],sparse(sum(A.A,2))))
    end
end

sum(A::Assoc) = sum(findnz(A.A)[3])

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
using PyCall
@pyimport matplotlib.pyplot as plt

function chart_weights(W, imsize; padding=0, annotation="", filename="", noshow=false, ordering=true)
    h, w = imsize
    n = size(W, 1)
    rows = round(Int,floor(sqrt(n)))
    cols = round(Int,ceil(n / rows))
    halfpad = div(padding, 2)
    dat = zeros(rows * (h + padding), cols * (w + padding))

    # Sort receptive fields by energy
    if ordering
        p = sum(W.^2,2)
        p = sortperm(vec(p);rev=true)
        W = W[p,:]
    end

    for i=1:n
        wt = W[i, :]
        wim = reshape(wt, imsize)
        r = div(i - 1, cols) + 1
        c = rem(i - 1, cols) + 1
        dat[(r-1)*(h+padding)+halfpad+1 : r*(h+padding)-halfpad,
            (c-1)*(w+padding)+halfpad+1 : c*(w+padding)-halfpad] = wim
    end

    normalize!(dat)

    return dat
end

function plot_hidden_activations(rbm::RBM,X::Mat{Float64})
    max_samples = 100
    n_samples = min(size(X,2),max_samples)
    x,_ = random_columns(X,n_samples)

    # Get all hidden activations for batch
    act = ProbHidCondOnVis(rbm,x)
    # Show this matrix of activations
    plt.imshow(act;interpolation="Nearest")
    plt.title("Hidden Unit Activations")
    plt.xlabel("Random Samples")
    plt.ylabel("Hidden Unit Index")
    plt.gray()

end

function plot_scores(mon::Monitor)
    ax_pl = plt.gca()
    ax_re = ax_pl[:twinx]()

    hpl = ax_pl[:plot](mon.Epochs,mon.PseudoLikelihood,"b^-",label="Pseudo-Likelihood")
    htl = ax_pl[:plot](mon.Epochs,mon.TAPLikelihood,"g^-",label="Tap-Likelihood")
    if mon.UseValidation
        hvpl = ax_pl[:plot](mon.Epochs,mon.ValidationPseudoLikelihood,"b^:",label="Pseudo-Likelihood (Validation)")
        hvtl = ax_pl[:plot](mon.Epochs,mon.ValidationTAPLikelihood,"g^:",label="Tap-Likelihood (Validation)")
    end
    ax_pl[:set_ylabel]("Normalized Likelihood")
    ax_pl[:set_ylim]((-0.3,0.0))

    hre = ax_re[:plot](mon.Epochs,mon.ReconError,"-*r",label="Recon. Error")
    if mon.UseValidation
        hvre = ax_re[:plot](mon.Epochs,mon.ValidationReconError,":*r",label="Recon. Error (Validation)")
    end
    ax_re[:set_ylabel]("Value")
    ax_re[:set_yscale]("log")

    plt.title("Scoring")
    plt.xlabel("Training Epoch")
    plt.xlim((1,mon.Epochs[mon.LastIndex]))
    plt.grid("on")
    if mon.UseValidation
        plt.legend(handles=[hpl;hvpl;htl;hvtl;hre;hvre],loc=3,fontsize=10)
    else
        plt.legend(handles=[hpl;htl;hre],loc=3,fontsize=10)
    end
end

function plot_evolution(mon::Monitor)
    hbt = plt.plot(mon.Epochs,mon.BatchTime_µs,"-k*",label="Norm. Batch time (µs)")

    plt.legend(handles=hbt,loc=1)
    plt.title("Evolution")
    plt.xlabel("Training Epoch")
    plt.xlim((1,mon.Epochs[mon.LastIndex]))
    plt.grid("on")
end

function plot_rf(rbm::RBM)
    # TODO: Implement RF display in the case of 1D signals
    rf = chart_weights(rbm.W,rbm.VisShape; padding=0,noshow=true)
    plt.imshow(rf;interpolation="Nearest")
    plt.title("Receptive Fields")
    plt.gray()
end

function plot_chain(rbm::RBM)
    # TODO: Implement Chain display in the case of 1D signals
    pc = chart_weights(rbm.persistent_chain_vis',rbm.VisShape; padding=0,noshow=true,ordering=false)
    plt.imshow(pc;interpolation="Nearest")
    plt.title("Visible Chain")
    plt.gray()
end

function plot_vbias(rbm::RBM)
    vectorMode = minimum(rbm.VisShape)==1 ? true : false

    if vectorMode
        plt.plot(rbm.vbias)
        plt.grid("on")
    else
        plt.imshow(reshape(rbm.vbias,rbm.VisShape);interpolation="Nearest")
    end
    plt.title("Visible Biasing")
    plt.gray()
end

function plot_weightdist(rbm::RBM)
    plt.hist(vec(rbm.W);bins=100)
    plt.title("Weight Distribution")
    plt.xlabel("Weight Value")
    plt.ylabel("Frequency")
end

function figure_refresh(figureHandle)
    figureHandle[:canvas][:draw]()
    plt.show(block=false)
    plt.pause(0.0001)
end

function WriteMonitorChartPDF(rbm::RBM,mon::Monitor,X::Mat{Float64},filename::AbstractString)
    savefig = plt.figure(5;figsize=(12,15))
    # Show Per-Epoch Progres
    savefig[:add_subplot](231)
        plot_scores(mon)

    savefig[:add_subplot](232)
        plot_evolution(mon)

    # Show Receptive fields
    savefig[:add_subplot](233)
        plot_rf(rbm)

    # Show the Visible chains/fantasy particle
    savefig[:add_subplot](234)
        plot_chain(rbm)

    # Show the current visible biasing
    savefig[:add_subplot](235)
        # plot_vbias(rbm)
        plot_hidden_activations(rbm,X)

    # Show the distribution of weight values
    savefig[:add_subplot](236)
        plot_weightdist(rbm)

    plt.savefig(filename;transparent=true,format="pdf",papertype="a4",frameon=true,dpi=300)
    plt.close()
end

function ShowMonitor(rbm::RBM,mon::Monitor,X::Mat{Float64},itr::Int;filename=[])
    fig = mon.FigureHandle

    if mon.MonitorVisual && itr%mon.MonitorEvery==0
        # Wipe out the figure
        fig[:clf]()

        # Show Per-Epoch Progres
        fig[:add_subplot](231)
            plot_scores(mon)

        fig[:add_subplot](232)
            plot_evolution(mon)

        # Show Receptive fields
        fig[:add_subplot](233)
            plot_rf(rbm)

        # Show the Visible chains/fantasy particle
        fig[:add_subplot](234)
            plot_chain(rbm)

        # Show the current visible biasing
        fig[:add_subplot](235)
            # plot_vbias(rbm)
            plot_hidden_activations(rbm,X)

        # Show the distribution of weight values
        fig[:add_subplot](236)
            plot_weightdist(rbm)

        figure_refresh(fig)
    end

    if mon.MonitorText && itr%mon.MonitorEvery==0
        li = mon.LastIndex
        ce = mon.Epochs[li]
        if mon.UseValidation
            @printf("[Epoch %04d] Train(pl : %0.3f, tl : %0.3f), Valid(pl : %0.3f, tl : %0.3f)  [%0.3f µsec/batch/unit]\n",ce,
                                                                                mon.TAPPseudoLikelihood[li],
                                                                                mon.TAPLikelihood[li],
                                                                                mon.ValidationPseudoLikelihood[li],
                                                                                mon.ValidationTAPLikelihood[li],
                                                                                mon.BatchTime_µs[li])
        else
            @printf("[Epoch %04d] Train(pl : %0.3f, tl : %0.3f)  [%0.3f µsec/batch]\n",ce,
                                                                           mon.TAPPseudoLikelihood[li],
                                                                           mon.TAPLikelihood[li],
                                                                           mon.BatchTime_µs[li])
        end
    end
end

# function chart_likelihood_evolution(pseudo, tap; filename="")

#     if length(filename) > 0
#         # Write to file if filename specified
#         LikelihoodPlot = plot(
#                                 layer(x=1:length(pseudo),y=pseudo,Geom.point, Geom.line),
#                                 layer(x=1:length(tap),y=tap, Geom.point, Geom.line, Theme(default_color=colorant"green")),
#                                 Guide.xlabel("epochs"),Guide.ylabel("Likelihood"),Guide.title("Evolution of likelihood for training set")
#                             )
#         draw(PDF(filename, 10inch, 6inch), LikelihoodPlot)
#     else
#         # Draw plot if no filename given
#        plot(
#                                 layer(x=1:length(pseudo),y=pseudo,Geom.point, Geom.line),
#                                 layer(x=1:length(tap),y=tap, Geom.point, Geom.line, Theme(default_color=colorant"green")),
#                                 Guide.xlabel("epochs"),Guide.ylabel("Likelihood"),Guide.title("Evolution of likelihood for training set")
#                             )
#     end
# end


text:
@testset "CoDaArray" begin
  table = (a=[1,2,3], b=[4,5,6], c=[7,8,9])
  array = CoDaArray(table)
  rows  = array
  @test length(array) == 3
  @test array[1] == Composition(a=1, b=4, c=7)
  @test array[2] == Composition(a=2, b=5, c=8)
  @test array[3] == Composition(a=3, b=6, c=9)
  @test Tables.istable(array)
  @test Tables.rowaccess(array) == true
  @test Tables.rows(array) == rows
  @test Tables.getcolumn(rows[1], 1) == 1.0
  @test Tables.getcolumn(rows[1], :a) == 1.0
  @test Tables.columnnames(rows[1]) == (:a, :b, :c)

  array = compose(jura, (:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn), keepcols=false)
  @test parts(array) == (:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn)
  @test length(array) == 359
  @test array[1] == Composition(Cd=1.74, Cu=25.72, Pb=77.36, Co=9.32, Cr=38.32, Ni=21.32, Zn=92.56)
  @test length(array.Cd) == 359
  @test length(array.Cu) == 359
  @test array.Cd[1] == 1.74
  @test array.Cu[1] == 25.72
  @test_throws KeyError array.INVALID

  table = compose(jura, (:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn))
  @test Tables.columnnames(table) == [:X,:Y,:Rock,:Land,:coda]
  @test Tables.getcolumn(table, :coda) == array

  table = compose(jura, (:Cd, :Cu, :Pb, :Co, :Cr, :Ni, :Zn), as=:comps)
  @test Tables.columnnames(table) == [:X,:Y,:Rock,:Land,:comps]
  @test Tables.getcolumn(table, :comps) == array
end


text:
function theme_cutie()
    p=plot(bg_color=:pink,fg_color = RGB(.7,.7,1),grid=(1, 0.9),size=(1e3,1e3),dpi=150,color = :red,markerstrokecolor = :red,key =false)
    return p
end

function theme_dark()
    p = Plots.plot(markerstrokewidth=.01,bg_color=:black,size=[1e3,1e3],dpi=150,xlabel = "x",
    ylabel = "y",zlabel = "z",fg_color = RGB(.7,.7,1),grid=(1, 0.9));
    return p
end
#aspect ratio need to be added later
theme_default = theme_dark


text:
using DashBootstrapComponents

email_input = dbc_formgroup([
    dbc_label("Email", html_for="example-email"),
    dbc_input(type="email", id="example-email", placeholder="Enter email"),
    dbc_formtext("Are you on email? You simply have to be these days", color="secondary"),
]);

password_input = dbc_formgroup([
    dbc_label("Password", html_for="example-password"),
    dbc_input(type="password", id="example-password", placeholder="Enter password"),
    dbc_formtext("A password stops mean people taking your stuff", color="secondary"),
]);

form = dbc_form([email_input, password_input]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


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

"""
    rc_filter(RC::Real)

Return a DSP.jl-compatible RC-filter.
"""
function rc_filter(RC::Real)
    T = float(typeof(RC))
    α = 1 / (1 + RC)
    Biquad(T(α), T(0), T(0), T(α - 1), T(0))
end
export rc_filter

"""
    cr_filter(RC::Real)

Return a DSP.jl-compatible CR-filter.
"""
function cr_filter(RC::Real)
    T = float(typeof(RC))
    α = RC / (RC + 1)
    Biquad(T(α), T(-α), T(0), T(-α), T(0))
end
export cr_filter

"""
    inv_cr_filter(RC::Real)

Return a DSP.jl-compatible inverse CR-filter.
"""
function inv_cr_filter(RC::Real)
    T = float(typeof(RC))
    α = 1 / (1 + RC)
    k = 1 + 1/RC
    Biquad(T(k), T(k * (α - 1)), T(0), T(-1), T(0))
end
export inv_cr_filter

"""
    crmod_filter(RC::Real)

Return a DSP.jl-compatible modified CR-filter.
"""
function crmod_filter(RC::Real)
    T = float(typeof(RC))
    α = RC / (RC + 1)
    Biquad(T(1), T(-1), T(0), T(-α), T(0))
end
export crmod_filter

"""
    inv_crmod_filter(RC::Real)

Return a DSP.jl-compatible inverse modified CR-filter.
"""
function inv_crmod_filter(RC::Real)
    T = float(typeof(RC))
    α = 1 / (1 + RC)
    Biquad(T(1), T(α - 1), T(0), T(-1), T(0))
end
export inv_crmod_filter

"""
    integrator_filter(gain::Real)

Return a DSP.jl-compatible integrator filter.
"""
function integrator_filter(gain::Real)
    T = float(typeof(gain))
    Biquad(T(gain), T(0), T(0), T(-1), T(0))
end
export integrator_filter

"""
    differentiator_filter(gain::Real)

Return a DSP.jl-compatible differentiator filter.
"""
function differentiator_filter(gain::Real)
    T = float(typeof(gain))
    Biquad(T(gain), T(-gain), T(0), T(0), T(0))
end
export differentiator_filter

"""
    integrator_cr_filter(gain::Real, RC::Real)

Return a DSP.jl-compatible integrator plus CR filter.
"""
function integrator_cr_filter(gain::Real, RC::Real)
    T = float(promote_type(typeof(gain), typeof(RC)))
    α = 1 / (1 + RC)
    Biquad(T(gain), T(-α), T(0), T(α - 1), T(0))
end
export integrator_cr_filter

"""
    integrator_crmod_filter(gain::Real, RC::Real)

Return a DSP.jl-compatible integrator plus modified CR filter.
"""
function integrator_crmod_filter(gain::Real, RC::Real)
    T = float(promote_type(typeof(gain), typeof(RC)))
    α = 1 / (1 + RC)
    Biquad(T(gain), T(0), T(0), T(α - 1), T(0))
end
export integrator_crmod_filter

"""
simple_csa_response_filter(τ_rise::Real, τ_decay::Real, gain::Real = one(τ_rise))

Return a DSP.jl-compatible filter that models the response of a typical
charge-sensitive amplifier (CSA).
"""
function simple_csa_response_filter(τ_rise::Real, τ_decay::Real, gain::Real = one(τ_rise))
    # TODO: Use a single biquad filter

    T = float(promote_type(promote_type(typeof(τ_rise), typeof(τ_decay)), typeof(gain)))
    rc_filter(T(τ_rise)) * integrator_cr_filter(T(gain), T(τ_decay))
end
export simple_csa_response_filter


text:
# WAHVectors
# ==========
#
# Bit Vectors compressed using a Word Aligned Hybrid method.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/WAHVectors.jl/blob/master/LICENSE.md

module WAHVectors

include("wah_elements.jl")
#include("31_bit_iteration.jl")
#include("vector_wahelement.jl")
#include("wah_vectors.jl")

#include("wah_cursor.jl")

export WAHVector

end


text:
using Peacock, PyPlot, LinearAlgebra
orthonormalise = Peacock.orthonormalise
get_k = Peacock.get_k
sample_path = Peacock.sample_path

a = 1; a1 = [1, 0]; a2 = [0, 1]
w = 0.17 # 柱子边长
d = w/2 # w/2 中间， a/2 - w/2 四周 a/4 中心

eps0 = 1; mu0 = 1
mu1 = [14 12.4im 0; -12.4im 14 0; 0 0 1]
eps1 = 15

function epf(x, y)
    if abs(x - d) < w/2 && abs(y - d) < w/2 || abs(x + d) < w/2 && abs(y - d) < w/2 || abs(x - d) < w/2 && abs(y + d) < w/2 || abs(x + d) < w/2 && abs(y + d) < w/2
        return eps1
    else
        return eps0
    end
end

function muf(x, y)
    if abs(x - d) < w/2 && abs(y - d) < w/2 || abs(x + d) < w/2 && abs(y - d) < w/2 || abs(x - d) < w/2 && abs(y + d) < w/2 || abs(x + d) < w/2 && abs(y + d) < w/2
        return mu1
    else
        return mu0
    end
end

#epf(x, y) = x^2 + y^2 < 0.2^2 ? 15 : 1
#muf(x, y) = x^2 + y^2 < 0.2^2 ? 1 : 1

d1 = d2 = 0.005
geometry = Geometry(epf, muf, a1, a2, d1, d2, TM)
solver = Peacock.FDFD.Solver(geometry, 2*[d1, d2])

kx = -pi:0.1:pi
ky = -pi:0.1:pi
modes = Peacock.FDFD.solve(solver,[0,0],TM,bands=1:1)
weighting = modes[1].weighting
data = zeros(ComplexF64, length(kx),length(ky)+1,2,length(weighting))
i = 1
for kx in kx
    j = 1
    for ky in ky
        for n = 1:2
            modes = Peacock.FDFD.solve(solver,[kx,ky],TM,bands=1:2)
            data0 = get_field_FDFD(modes[n].data, modes[n].basis,k0=modes[n].k0)
            data0 = data0[:] / sqrt(abs(dot(data0[:], weighting.*data0[:])))
            data[i,j,n,:] = data0[:]
        end
        j = j + 1
    end
    data[i,end,1,:] = data[i,1,1,:]
    data[i,end,2,:] = data[i,1,2,:]
    i = i + 1
end

vy = zeros(ComplexF64, length(kx),length(ky),2)
W = zeros(ComplexF64, 2,2)
states = zeros(ComplexF64, length(kx),length(ky),2,2)
for i = 1:length(kx)
    for j1 = 1:length(ky)
        F = I
        for j2 = 0:length(ky)-1
            if j1 + j2 <= length(ky)
                j = j1+j2
            else
                j = j1 + j2 - length(ky)
            end
            for k = 1:2
                for l = 1:2
                    W[k,l] = dot(data[i,j,k,:],weighting.*data[i,j+1,l,:])
                end
            end
            F = F * W
        end
        vals,state = eigen(F,sortby=angle)
        vy[i,j1,:] = angle.(vals)
        states[i,j1,:,:] = state
    end
end

wx = zeros(ComplexF64, length(kx)+1,length(ky),length(weighting))
for j = 1:length(ky) #ky
    for i = 1:length(kx) #kx
        wx[i,j,:] = states[i,j,1,1] * data[i,j,1,:] + states[i,j,2,1] * data[i,j,2,:]
        wx[i,j,:] = wx[i,j,:]/sqrt(abs(dot(wx[i,j,:],weighting.*wx[i,j,:])))
    end
    wx[end,j,:] = wx[1,j,:]
end

px = zeros(length(ky))
for j = 1:length(ky) #ky
    W = 1
    for i = 1:length(kx) #kx
        W = W*dot(wx[i,j,:],weighting.*wx[i+1,j,:])
    end
#    vals = eigvals(W)
    px[j] = angle.(W)
#    px[j] = imag(log(W/2/pi))
end
plot(px,"*")
ylim(-pi,pi)
#figure()
#plot(vy,"*")


text:
# change the loop order in Base.
function naive_mul!(o::AbstractMatrix{T0}, a::AbstractMatrix{T1}, b::AbstractMatrix{T2}, α=one(T0), β=zero(T0)) where {T0,T1,T2}
    @assert size(a, 2) == size(b, 1) && size(o) == (size(a, 1), size(b, 2))
    a = convert(Matrix, a)
    b = convert(Matrix, b)
    for j=1:size(b, 2)
        if !iszero(β)
            @inbounds for i=1:size(a, 1)
                o[i,j] = β * o[i,j]
            end
        else
            @inbounds for i=1:size(a, 1)
                o[i,j] = zero(T0)
            end
        end
        for k=1:size(a, 2)
            for i=1:size(a, 1)
                @inbounds o[i,j] += α * a[i,k] * b[k,j]
            end
        end
    end
    return o
end

# Overwrite the `mul!` in LinearAlgebra (also changes the behavior of `*` in Base)!
for TA in [:(AbstractMatrix{T} where T<:TropicalTypes), :(Transpose{T,S} where {T<:TropicalTypes,S<:AbstractVecOrMat{T}})]
    for TB in [:(AbstractMatrix{T} where T<:TropicalTypes), :(Transpose{T,S} where {T<:TropicalTypes,S<:AbstractVecOrMat{T}})]
        @eval @inline function LinearAlgebra.mul!(o::AbstractMatrix{TO}, a::$TA, b::$TB, α::Number, β::Number) where TO
            α = _convert_to_static(TO, α)
            β = _convert_to_static(TO, β)
            naive_mul!(o, a, b, α, β)
        end
    end
end

Base.:*(a::T, b::StaticInt{0}) where T<:TropicalTypes = zero(T)
Base.:*(a::T, b::StaticInt{1}) where T<:TropicalTypes = a
Base.:*(b::StaticInt{0}, a::T) where T<:TropicalTypes = zero(T)
Base.:*(b::StaticInt{1}, a::T) where T<:TropicalTypes = a


text:
import Pkg
Pkg.add("PyCall")
Pkg.add("BitInformation")
Pkg.add("StatsBase")
Pkg.status()


text:
# function create_description()
# end

descript(::Any, i) = nothing
descript(::Type{T}, i) where {T <: AbstractParameter} = (T,LengthParamDescription(T, i))

function type_to_curlycall(::Type{T}) where {T}
    q = Expr(:curly, Symbol(T.name))
    for i ∈ eachindex(T.parameters)
        push!(q.args, T.parameters[i])
    end
    Expr(:call, q)
end

function deconstruct_namedtuple(::Type{NamedTuple{N,T}}) where {N,T}
    names = Symbol[N...]
    Tp = T.parameters
    lengthdescripts = LengthParamDescription[]
    descripts = Expr[]
    track = Vector{Bool}(undef, length(names))
    j = 1
    for i ∈ eachindex(names)
        dld = descript(Tp[i], j)
        track[i] = notnothing = !isnothing(dld)
        if notnothing
            push!(descripts, type_to_curlycall(first(dld)))
            push!(lengthdescripts, last(dld))
            j += 1
        end
    end
    names, track, lengthdescripts, descripts
end

preprocess(m::Model, descriptnt) = preprocess!(copy(m), descriptnt)
function preprocess!(m::Model, descriptnt)
    names, track, lengthdescripts, descripts = deconstruct_namedtuple(descriptnt)
    preprocess!(m, names, track, lengthdescripts, descripts)
end

function preprocess!(m::Model, names::Vector{Symbol}, track::Vector{Bool}, lengthdescripts::Vector{LengthParamDescription}, descripts)
    ReverseDiffExpressions.reset_var_tracked!(m)
    offsets = DistributionParameters.parameter_offsets(lengthdescripts)
    parametervar = addvar!(m, Symbol("#θ#")); parametervar.tracked = true
    datavar = addvar!(m, Symbol("#DATA#")); parametervar.tracked = false
    target = m.vars[0]
    j = 0
    for i ∈ eachindex(names)
        name = names[i]
        var = getvar!(m, name)
        if (var.tracked = track[i]) # initialized by loading
            descriptvar = addvar!(m, Symbol(""))
            offsetvar = addvar!(m, Symbol(""))
            constrainret = addvar!(m, gensym(:targetconstrained))
            descriptvar.ref = descripts[(j += 1)]
            offsetvar.ref = offsets[j]
            constrainfunc = Func(Instruction(:DistributionParameters,:constrain), false)
            uses!(constrainfunc, parametervar)
            uses!(constrainfunc, offsetvar)
            uses!(constrainfunc, descriptvar)
            returns!(constrainfunc, constrainret)
            addfunc!(m, constrainfunc)
            getconstindex!(m, target, constrainret, 1)
            getconstindex!(m, var, constrainret, 2)
        else
            getprop = Func(Instruction(:Base,:getproperty), false)
            namevar = addvar!(m, Symbol("")); namevar.ref = name
            uses!(getprop, datavar)
            uses!(getprop, namevar)
            returns!(getprop, var)
            addfunc!(m, getprop)
        end
    end
    parametervar.tracked = true
    ReverseDiffExpressions.propagate_var_tracked!(m)
    m
end


text:
function readTestInput()
  draw = [7,4,9,5,11,17,23,2,0,14,21,24,10,16,13,6,15,25,12,22,18,20,8,19,3,26,1]
  boards = [
    [22 13 17 11  0;
    8  2 23  4 24;
    21  9 14 16  7;
    6 10  3 18  5;
    1 12 20 15 19],

    [3 15  0  2 22;
    9 18 13 17  5;
    19  8  7 25 23;
    20 11 10 24  4;
    14 21 16 12  6],

    [14 21 17 24  4;
    10 16 15  9 19;
    18  8 23 26 20;
    22 11 13  6  5;
    2  0 12  3  7]
  ]
  return draw,boards
end

function readInput()
  open(joinpath(@__DIR__, "./4.txt")) do io
    # read draw
    draw = parse.(Int,split(strip(readline(io)),","))

    # ignore the next blank line
    readline(io)

    # read the boards into 5x5 matrices
    data = split.(split(read(io,String), "\n\n"),"\n")
    data = map(x->split.(x," "),data)
    data = map(x->map(y->filter(z->z!="",y),x),data)
    boards = map(x->parse.(Int,hcat(x...)),data)
    return draw,boards
  end
end

function score(x,y)
  if x == -1 && y == -1
    return 0
  elseif x == -1
    return y
  elseif y == -1
    return x
  else
    return x+y
  end
end

function part1(draw,boards)
  for num in draw, board in boards
    board[board.==num] .= -1
    if sum(all(i->i==-1,board,dims=1)) > 0 || sum(all(i->i==-1,board,dims=2)) > 0
      return board, num
    end
  end
end

function part2(draw,boards)
  winning_board = undef
  winning_num = undef
  while !isempty(boards)
    winning_board, winning_num = part1(draw,boards)
    setdiff!(boards,[winning_board])
  end
  return winning_board, winning_num
end

draw,boards = readInput()
winning_board, winning_num = part1(draw,boards)
board_score = reduce(score, winning_board)
board_score * winning_num

draw,boards = readInput()
winning_board, winning_num = part2(draw,boards)
board_score = reduce(score, winning_board)
board_score * winning_num


text:
module TestTraits

# using Revise
using MLJBase
import MLJBase
using Test

mutable struct Dummy <: Deterministic{Int} end

MLJBase.target_kind(::Type{Dummy}) = :numeric
MLJBase.target_quantity(::Type{Dummy}) = :univariate
MLJBase.inputs_can_be(::Type{Dummy}) = [:numeric, ]
MLJBase.is_pure_julia(::Type{Dummy}) = :yes
MLJBase.package_name(::Type{Dummy}) = "MultivariateStats"
MLJBase.package_uuid(::Type{Dummy}) = "6f286f6a-111f-5878-ab1e-185364afe411"
MLJBase.package_url(::Type{Dummy}) = "https://mickey.mouse.org"

d = Dict(:model_name => "TestTraits.Dummy",
         :is_pure_julia => :yes,
         :package_uuid  => "6f286f6a-111f-5878-ab1e-185364afe411",
         :package_name  => "MultivariateStats",
         :target_is     => Symbol[:deterministic, :numeric, :univariate],
         :inputs_can_be => Symbol[:numeric],
         :package_url   => "https://mickey.mouse.org")

@show info(Dummy)[:model_name]
# @test info(Dummy) == d
# @test info(Dummy()) == d

end
true


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiMetricsV1beta1PodMetricsList <: SwaggerModel
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ListMeta } # spec name: metadata
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiMetricsV1beta1PodMetrics} } # spec name: items

    function IoK8sApiMetricsV1beta1PodMetricsList(;metadata=nothing, items=nothing)
        o = new()
        validate_property(IoK8sApiMetricsV1beta1PodMetricsList, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiMetricsV1beta1PodMetricsList, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        o
    end
end # type IoK8sApiMetricsV1beta1PodMetricsList

const _property_map_IoK8sApiMetricsV1beta1PodMetricsList = Dict{Symbol,Symbol}(Symbol("metadata")=>Symbol("metadata"), Symbol("items")=>Symbol("items"))
const _property_types_IoK8sApiMetricsV1beta1PodMetricsList = Dict{Symbol,String}(Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ListMeta", Symbol("items")=>"Vector{IoK8sApiMetricsV1beta1PodMetrics}")
Base.propertynames(::Type{ IoK8sApiMetricsV1beta1PodMetricsList }) = collect(keys(_property_map_IoK8sApiMetricsV1beta1PodMetricsList))
Swagger.property_type(::Type{ IoK8sApiMetricsV1beta1PodMetricsList }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiMetricsV1beta1PodMetricsList[name]))}
Swagger.field_name(::Type{ IoK8sApiMetricsV1beta1PodMetricsList }, property_name::Symbol) =  _property_map_IoK8sApiMetricsV1beta1PodMetricsList[property_name]

function check_required(o::IoK8sApiMetricsV1beta1PodMetricsList)
    true
end

function validate_property(::Type{ IoK8sApiMetricsV1beta1PodMetricsList }, name::Symbol, val)
end


text:
########################################################################
function ppca(X; Q=1)
########################################################################

    D, N = size(X)

    μ = vec(mean(X, dims=2))

    JITTER = 1e-8

    @printf("Running PPCA for %d number of data items of dimension %d, projected to %d components\n", N, D, Q)

    #---------------------------------------------
    function marginalLogLikelihood(W, σ)
    #---------------------------------------------

        sum(logpdf(MvNormal(μ, W*W' + σ*σ*I + JITTER*I), X))

    end

    #---------------------------------------------
    function unpack(param)
    #---------------------------------------------

        local W = reshape(param[1:end-1], D, Q)

        local σ = exp(param[end])

        return W, σ

    end

    #---------------------------------------------
    function objective(param)
    #---------------------------------------------

        @assert(length(param) == D*Q + 1)

        local W, σ = unpack(param)

        return -1.0 * marginalLogLikelihood(W, σ)

    end

    #---------------------------------------------
    # Run optimiser
    #---------------------------------------------

    opt    = Optim.Options(show_trace = true, iterations = 10_000)

    result = optimize(objective, [randn(D*Q); randn()*3], LBFGS(), opt, autodiff=:forward)

    W, σ   = unpack(result.minimizer)

    @show W, μ, σ

    #---------------------------------------------
    # Define projections
    #---------------------------------------------

    M = W'*W  + σ*σ*I # eq. (12.41) in Bishop

    proj(x::Array{T,1} where T<:Real) = vec(M \ (W'*(x - μ)))[1]

    fwd(z::Real) = vec(W*z + μ)

    return W, μ, σ, proj, fwd

end


text:
function monomial1d(x::Array{Float64,2},p::Int)::Tuple{Array{Float64,2},Array{Float64,2}}

n = length(x);
v0 = zeros(Float64,n,1);
v1 = ones(Float64,n,1);

if (p==0)
    f = v1;
    fx = v0;
elseif (p==1)
    f = [v1 x];
    fx = [v0 v1];
elseif (p==2)
    f = [v1 x x.^2];
    fx = [v0 v1 2*x];
elseif (p==3)
    f = [v1 x x.^2 x.^3];
    fx = [v0 v1 2*x 3*x.^2];
elseif (p==4)
    f = [v1 x x.^2 x.^3 x.^4];
    fx = [v0 v1 2*x 3*x.^2 4*x.^3];
else
    error("polynomial degree must be less than 5");
end

return f, fx;

end

function monomial2d(xy::Array{Float64,2},p::Int)::Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2}}

n = size(xy,1);
v0 = zeros(Float64,n,1);
v1 = ones(Float64,n,1);

x = xy[:,1];
y = xy[:,2];

if (p==0)
    f = v1;
    fx = v0;
    fy = v0;
elseif (p==1)
    f = [v1 x y];
    fx = [v0 v1 v0];
    fy = [v0 v0 v1];
elseif (p==2)
    f = [v1 x y x.*y x.^2 y.^2];
    fx = [v0 v1 v0 y 2*x v0];
    fy = [v0 v0 v1 x v0 2*y];
elseif (p==3)
    f = [v1 x y x.*y x.^2 y.^2 x.*y.^2 (x.^2).*y x.^3 y.^3];
    fx = [v0 v1 v0 y 2*x v0 y.^2 2*x.*y 3*x.^2 v0];
    fy = [v0 v0 v1 x v0 2*y 2*x.*y x.^2 v0 3*y.^2];
elseif (p==4)
    f = [v1 x y x.*y x.^2 y.^2 x.*y.^2 (x.^2).*y x.^3 y.^3 x.*y.^3 (x.^2).*y.^2 (x.^3).*y x.^4 y.^4];
    fx = [v0 v1 v0 y 2*x v0 y.^2 2*x.*y 3*x.^2 v0 y.^3 2*x.*y.^2 3*(x.^2).*y 4*x.^3 v0];
    fy = [v0 v0 v1 x v0 2*y 2*x.*y x.^2 v0 3*y.^2 3*x.*y.^2 2*(x.^2).*y x.^3 v0 4*y.^3];
else
    error("polynomial degree must be less than 5");
end

return f, fx, fy;

end

function monomial3d(xyz::Array{Float64,2},p::Int)::Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}

n = size(xyz,1);
v0 = zeros(Float64,n,1);
v1 = ones(Float64,n,1);

x = xyz[:,1];
y = xyz[:,2];
z = xyz[:,3];

if (p==0)
    f = v1;
    fx = v0;
    fy = v0;
    fz = v0;
elseif (p==1)
    f = [v1 x y z];
    fx = [v0 v1 v0 v0];
    fy = [v0 v0 v1 v0];
    fz = [v0 v0 v0 v1];
elseif (p==2)
    f = [v1 x y z x.*y x.*z y.*z x.^2 y.^2 z.^2];
    fx = [v0 v1 v0 v0 y z v0 2*x v0 v0];
    fy = [v0 v0 v1 v0 x v0 z v0 2*y v0];
    fz = [v0 v0 v0 v1 v0 x y v0 v0 2*z];
elseif (p==3)
    f = [v1 x y z x.*y x.*z y.*z x.^2 y.^2 z.^2 x.*y.*z x.*y.^2 x.*z.^2 (x.^2).*y (x.^2).*z y.*z.^2 (y.^2).*z x.^3 y.^3 z.^3];
    fx = [v0 v1 v0 v0 y z v0 2*x v0 v0 y.*z y.^2 z.^2 2*x.*y 2*x.*z v0 v0 3*x.^2 v0 v0];
    fy = [v0 v0 v1 v0 x v0 z v0 2*y v0 x.*z 2*x.*y v0 x.^2 v0 z.^2 2*y.*z v0 3*y.^2 v0];
    fz = [v0 v0 v0 v1 v0 x y v0 v0 2*z x.*y v0 2*x.*z v0 x.^2 2*y.*z y.^2 v0 v0 3*z.^2];
elseif (p==4)
    f = [v1 x y z x.*y x.*z y.*z x.^2 y.^2 z.^2 ...
         x.*y.*z x.*y.^2 x.*z.^2 (x.^2).*y (x.^2).*z y.*z.^2 (y.^2).*z x.^3 y.^3 z.^3 ...
         x.*y.*z.^2 (x.*y.^2).*z (x.^2).*y.*z (x.^2).*y.^2 (x.^2).*z.^2 (y.^2).*z.^2 ...
         (x.^3).*y (x.^3).*z x.*y.^3 (y.^3).*z x.*z.^3 y.*z.^3 x.^4 y.^4 z.^4];
    fx = [v0 v1 v0 v0 y z v0 2*x v0 v0 y.*z y.^2 z.^2 2*x.*y 2*x.*z v0 v0 3*x.^2 v0 v0...
          y.*z.^2 (y.^2).*z 2*x.*y.*z 2*x.*y.^2 2*x.*z.^2 v0 ...
          3*(x.^2).*y 3*(x.^2).*z y.^3 v0 z.^3 v0 4*x.^3 v0 v0];
    fy = [v0 v0 v1 v0 x v0 z v0 2*y v0 x.*z 2*x.*y v0 x.^2 v0 z.^2 2*y.*z v0 3*y.^2 v0...
          x.*z.^2 2*x.*y.*z (x.^2).*z 2*(x.^2).*y v0 2*y.*z.^2 ...
          x.^3 v0 3*x.*y.^2 3*(y.^2).*z v0 z.^3 v0 4*y.^3 v0];
    fz = [v0 v0 v0 v1 v0 x y v0 v0 2*z x.*y v0 2*x.*z v0 x.^2 2*y.*z y.^2 v0 v0 3*z.^2 ...
          2*x.*y.*z x.*y.^2 (x.^2).*y v0 2*(x.^2).*z 2*(y.^2).*z ...
          v0 x.^3 v0 y.^3 3*x.*z.^2 3*y.*z.^2 v0 v0 4*z.^3];
else
    error("polynomial degree must be less than 5");
end

return f, fx, fy, fz;

end

function simplexmonomial(x::Array{Float64,2},porder::Int)::Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}

dim=size(x,2);

if dim==1
     f,fx=monomial1d(x,porder);      # 1D
     fy=[0.0 0.0];
     fz=[0.0 0.0];
elseif dim==2
     f,fx,fy=monomial2d(x,porder);   # 2D
     fz=[0.0 0.0];
elseif dim==3
     f,fx,fy,fz=monomial3d(x,porder);# 3D
else
     error("Only can handle dim=1, dim=2 or dim=3");
end

return f, fx, fy, fz;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# SteinDistribution
# An abstract type representing a distribution P on R^d with support
# (alpha_1,beta_1) x ... x (alpha_d,beta_d) with -inf <= alpha_j < beta_j <= inf
abstract SteinDistribution
import Base.rand

# Returns first Stein factor
function getC1(d::SteinDistribution)
    return d.c1
end
# Returns second Stein factor
function getC2(d::SteinDistribution)
    return d.c2
end
# Returns third Stein factor
function getC3(d::SteinDistribution)
    return d.c3
end
# Returns lower bound of support for j-th dimension marginal of P
# Should be overridden for constrained distributions
function supportlowerbound(d::SteinDistribution, j::Int64)
    -Inf
end
# Returns upper bound of support for j-th dimension marginal of P
# Should be overridden for constrained distributions
function supportupperbound(d::SteinDistribution, j::Int64)
    Inf
end

# a convenience method if this is only defined for 1d vectors
function gradlogdensity(dist::SteinDistribution,
                        xs::Array{Float64,2})
    (n, d) = size(xs)
    if n < 1
        return []
    end
    gradlogdensities = Array(Float64, n, d)
    for ii in 1:n
        gradlogdensities[ii,:] = gradlogdensity(dist, vec(xs[ii,:]))
    end
    gradlogdensities
end


text:
using VarInfoX
using BenchmarkTools

SUITE = BenchmarkGroup()

SUITE["seq"] = @benchmarkable VarInfoX.varinfo_seq(Base; recursive = true, all = true)
# SUITE["pr42123"] =
#     @benchmarkable VarInfoX.varinfo_pr42123(Base; recursive = true, all = true)
SUITE["dac_names"] =
    @benchmarkable VarInfoX.varinfo_dac_names(Base; recursive = true, all = true)
SUITE["dac_module"] =
    @benchmarkable VarInfoX.varinfo_dac_module(Base; recursive = true, all = true)
SUITE["parallel_names"] =
    @benchmarkable VarInfoX.varinfo_parallel_names(Base; recursive = true, all = true)


text:
using LinearAlgebra:
        Diagonal,
        copyto!,
        diag,
        I

struct CoupledResonators
    yr::AbstractMatrix
    yi::AbstractMatrix
    wr::AbstractMatrix
    wi::AbstractMatrix
    γin::AbstractMatrix
    γex::AbstractMatrix
end

function CoupledResonators()
    return CoupledResonators(
        zeros(0, 0),
        zeros(0, 0),
        zeros(0, 0),
        zeros(0, 0),
        Diagonal{Real}([]),
        Diagonal{Real}([]),
     )
end

modes(cr::CoupledResonators) = size(cr.yr, 1)

function addModes(cr::CoupledResonators, l::Int)
    if l == 0
        return cr
    end
    return CoupledResonators(
    Matrix(cr.yr ⊕ zeros(l, l)),
    Matrix(cr.yi ⊕ zeros(l, l)),
    Matrix(cr.wr ⊕ zeros(l, l)),
    Matrix(cr.wi ⊕ zeros(l, l)),
    Diagonal( [ diag(cr.γin)..., zeros(l)... ] ),
    Diagonal( [ diag(cr.γex)..., ones(l)... ] )
    )
end

# add g a⁺ⱼaₖ + conj(g) aⱼa⁺ₖ
function addPassive(cr::CoupledResonators,
    g::Number, j::Int, k::Int)
    m = size(cr.yr, 1)
    cr = addModes(cr, max(j - m, k - m, 0 ))
    yr = cr.yr
    yi = cr.yi
    wi = cr.wi
    wr = cr.wr
    γin = cr.γin
    γex = cr.γex
    yr[j, k] = real(g)
    yr[k, j] = real(g)
    yi[j, k] = imag(g)
    yi[k, j] = -imag(g)
    return CoupledResonators(yr, yi, wr, wi, γin, γex)
end

# add g a⁺ⱼa⁺ₖ + conj(g) aⱼaₖ
function addActive(cr::CoupledResonators,
    g::Number, j::Int, k::Int)
    m = size(cr.yr, 1)
    cr = addModes(cr, max(j - m, k - m, 0 ))
    yr = cr.yr
    yi = cr.yi
    wi = cr.wi
    wr = cr.wr
    γin = cr.γin
    γex = cr.γex
    wr[j, k] = real(g) / 2
    wr[k, j] = real(g) / 2
    wi[j, k] = imag(g) / 2
    wi[k, j] = imag(g) / 2
    return CoupledResonators(yr, yi, wr, wi, γin, γex)
end

function addGammaIn(cr::CoupledResonators,
    γ::Number, l::Int)
    m = size(cr.yr, 1)
    cr = addModes(cr, max(l - m, 0 ))
    yr = cr.yr
    yi = cr.yi
    wi = cr.wi
    wr = cr.wr
    γin = cr.γin
    γex = cr.γex
    γin[l, l] = γ
    return CoupledResonators(yr, yi, wr, wi, γin, γex)
end

function addGammaEx(cr::CoupledResonators,
    γ::Number, l::Int)
    m = size(cr.yr, 1)
    cr = addModes(cr, max(l - m, 0 ))
    yr = cr.yr
    yi = cr.yi
    wi = cr.wi
    wr = cr.wr
    γin = cr.γin
    γex = cr.γex
    γex[l, l] = γ
    return CoupledResonators(yr, yi, wr, wi, γin, γex)
end

const e0 = I(4)

const e1 = [0 0 0 1;
            0 0 1 0;
            0 1 0 0;
            1 0 0 0]

const e2 = [0 0 1 0;
            0 0 0 -1;
            1 0 0 0;
            0 -1 0 0]

const e3 = [1 0 0 0;
            0 1 0 0;
            0 0 -1 0;
            0 0 0 -1]

function scatteringMatrix(ω, yr, yi, wr, wi, γin, γex)
    m = size(yr, 1)
    Θ = ω * inv(γex + γin)
    Yr = inv(sqrt.(γex)) * yr * inv(sqrt.(γex))
    Yi = inv(sqrt.(γex)) * yi * inv(sqrt.(γex))
    Wr = inv(sqrt.(γex)) * wr * inv(sqrt.(γex))
    Wi = inv(sqrt.(γex)) * wi * inv(sqrt.(γex))
    Γ = inv(sqrt.(γex)) * (γex + γin) * inv(sqrt.(γex)) / 2
    S = e0 ⊗ I(m) - inv(
          e0 ⊗ (Yi + Γ) -
          e1 ⊗ Wr +
          e2 ⊗ Wi -
          (e1 * e2) ⊗ Yr +
          (e1 * e2 * e3) ⊗ Θ
    )
    order = vcat([ [i, i + m] for i in 1:m]...,  [ [i + 2*m, i + 3*m] for i in 1:m]...)
    return S[order, order]
end
function scatteringMatrix(ω::Number, cr::CoupledResonators)
    yr = cr.yr
    yi = cr.yi
    wr = cr.wr
    wi = cr.wi
    γin = cr.γin
    γex = cr.γex
    return scatteringMatrix(ω, yr, yi, wr, wi, γin, γex)
end


text:
# UNDER CONSTRUCTION

mutable struct Parallel_row
    ind::Int64
    ratio::Float64
    u::Float64
    g::Float64

    function Parallel_row(ind::Int64,ratio::Float64,u::Float64,g::Float64)
        return new(ind,ratio,u,g)
    end
end

mutable struct Parallel_row_group
    ls::Array{Parallel_row,1}
    first::Int64
    u::Float64 # combined u
end

# find and detect all parrell constraints
# e.g., a' * x and θ * a' * x
# not only eliminate duplicate bound constraints
# algorithm: scale all vectors so first element is one and sort list
type Clever_Symmetric_KKT_solver <: abstract_KKT_system_solver
    # abstract_KKT_system_solver
    ls_solver::abstract_linear_system_solver
    factor_it::Class_iterate
    delta_x_vec::Array{Float64,1}
    delta_s_vec::Array{Float64,1}
    rhs::System_rhs
    dir::Class_point
    kkt_err_norm::Class_kkt_error
    rhs_norm::Float64
    pars::Class_parameters
    schur_diag::Array{Float64,1}
    true_x_diag::Array{Float64,1}

    ready::Symbol

    # Symmetric_KKT_solver only
    Q::SparseMatrixCSC{Float64,Int64}

    # these are used to eliminate parrellel constraints from the linear system solve
    para_row_info::Array{Parallel_row_group,1}
    first_para_indicies::Array{Int64,1}

    function Clever_Symmetric_KKT_solver()
      return new()
    end
end

function initialize!(kkt_solver::Clever_Symmetric_KKT_solver, intial_it::Class_iterate)
    # call this before running using the kkt_solver
    initialize!(kkt_solver.ls_solver)
    kkt_solver.dir = zero_point(dim(intial_it),ncon(intial_it))

    #start_advanced_timer(timer, "symmetric/form_system/compute_indicies");
    kkt_solver.first_para_indicies, kkt_solver.para_row_info = compute_indicies(get_jac(intial_it))
    #pause_advanced_timer(timer, "symmetric/form_system/compute_indicies");
end

function columns_are_same(A::SparseMatrixCSC{Float64,Int64},i::Int64,j::Int64)
    a_i = A[:,i]
    a_j = A[:,j]
    if a_i.nzind == a_j.nzind
        if length(a_i.nzval) == length(a_j.nzval)
            if length(a_i.nzval) > 0
                ratio = a_i.nzval[1] / a_j.nzval[1]
                if norm(a_i - a_j * ratio,2) < 1e-16
                    return true
                else
                    return false
                end
            else
                warn("zero elements in column")
                return true
            end
        else
            return false
        end
    else
        return false
    end
end

function rescale_cols(A::SparseMatrixCSC{Float64,Int64})
    A_rescaled = deepcopy(A)
    new_vals = A_rescaled.nzval
    i = 1
    for col = 1:size(A,2)
        a_col = A[:,col]
        if length(a_col.nzind) > 0
            tmp = a_col.nzval / a_col.nzval[1]
            l = length(tmp)
            new_vals[i:(i+l-1)] = tmp
            i += l
        end
    end

    return A_rescaled
end

function compare_columns(A::SparseMatrixCSC{Float64,Int64},i::Int64,j::Int64)
    const a_i_nzind = A[:,i].nzind
    const a_j_nzind = A[:,j].nzind
    if length(a_j_nzind) == 0
        return false
    elseif length(a_i_nzind) == 0
        return true
    end

    a_i_first_ind = a_i_nzind[1]
    a_j_first_ind = a_j_nzind[1]

    if a_i_first_ind != a_j_first_ind
        if a_j_first_ind < a_i_first_ind
            return true
        else
            return false
        end
    else
        if length(a_i_nzind) < length(a_j_nzind)
            return true
        elseif length(a_i_nzind) > length(a_j_nzind)
            return false
        end

        for w = 1:length(a_i_nzind)
            if a_i_nzind[w] < a_j_nzind[w]
                return true
            elseif a_i_nzind[w] > a_j_nzind[w]
                return false
            end
        end

        a_i = A[:,i]
        a_j = A[:,j]
        for k in a_i.nzind
            if a_i[k] < a_j[k]
                return true
            elseif a_i[k] > a_j[k]
                return false
            end
        end
        if i < j
            return true
        else
            return false
        end
    end
end

function sorted_col_list(A::SparseMatrixCSC{Float64,Int64})
    m = size(A,2)
    rescaled_A = rescale_cols(A)
    function comp(i::Int64,j::Int64)
        return compare_columns(rescaled_A,i,j)
    end
    sorted_cols = collect(1:m)
    sort!(sorted_cols,lt=comp)
    return sorted_cols
end

function compute_breakpoints(A::SparseMatrixCSC{Float64,Int64},sorted_cols::Array{Int64,1})
    break_points = Array{Int64,1}()
    for bp = 1:length(sorted_cols)
        if bp == 1
            push!(break_points,bp)
        else
            prev = sorted_cols[bp-1]
            cur = sorted_cols[bp]
            if !columns_are_same(A,prev,cur)
                push!(break_points,bp)
            end
        end
    end
    return break_points
end

function compute_indicies(J::SparseMatrixCSC{Float64,Int64})
    dic = Dict{Any,Array{Int64,1}}()
    J_T = J'
    # REWRITE FASTER BY EMPLOYING SORTING FUNCTION.
    # a < b if
    # ahat = a / a[a_first_index], bhat = b / b[b_first_index]
    # find first i with ahat[i] < bhat[i]
    sorted_cols = sorted_col_list(J_T)

    break_points = compute_breakpoints(J_T,sorted_cols)

    no_para_indicies = sort(sorted_cols[break_points])
    para_row_info = Array{Parallel_row_group,1}()
    for k = 1:length(break_points)
        bp = break_points[k]
        if k < length(break_points)
            end_bp = break_points[k+1]
        else
            end_bp = length(sorted_cols)+1
        end
        first = sorted_cols[bp]
        ind_ls = sorted_cols[bp:(end_bp-1)]
        ls = Array{Parallel_row,1}()
        for i in ind_ls
            if length(J_T[:,i].nzval) > 0
                ratio = J_T[:,i].nzval[1] / J_T[:,ind_ls[1]].nzval[1]
            else
                ratio = 1.0
            end
            push!(ls, Parallel_row(i,ratio,NaN,NaN))
        end

        prg = Parallel_row_group(ls,first,NaN)
        push!(para_row_info, prg)
    end

    sort!(para_row_info,by=x->x.first)
    return no_para_indicies, para_row_info
end

function check_para_rows()
    # TODO
end

function compute_indicies(J::SparseMatrixCSC{Float64,Int64}, diag_vals::Array{Float64,1})
    ##################################################################################################################
    # INPUT
    # J::SparseMatrixCSC{Float64,Int64}
    # diag_vals::Array{Float64,1}
    #
    # OUTPUT
    #  no_para_indicies::Array{Int64,1} is a sorted list of rows that are not parallel to any other rows
    #  para_row_info::Array{} is a list containing each group of rows that are parallel.
    ##################################################################################################################

    no_para_indicies, para_row_info = compute_indicies(J)
    update_indicies!(para_row_info,diag_vals)
    return no_para_indicies, para_row_info
end

function update_indicies!(para_row_info::Array{Parallel_row_group,1}, diag_vals::Array{Float64,1})
    for row_group in para_row_info
        ls = row_group.ls
        u_inv = 0.0
        for row in ls
            row.u = diag_vals[row.ind]
            u_inv += (row.ratio)^2 * (row.u)^(-1.0)
        end

        group_u = 1.0/u_inv
        if (group_u <= 0.0)
            error("$group_u == group_u <= 0.0")
        end
        row_group.u = group_u

        for row in ls
            row.g = group_u * row.ratio * (row.u)^(-1.0)
        end
    end
end

#=
A = J[lb,:]
B = J[ineq,:]
D_l, D_u, D_eq

D_bd = (D_1.^(-1) + D_2.^(-1))^(-1)

H A' B'
A D_bd 0
B 0 D_eq
=#

function form_system!(kkt_solver::Clever_Symmetric_KKT_solver, iter::Class_iterate, timer::class_advanced_timer)
    k = kkt_solver
    start_advanced_timer(timer, "symmetric/form_system");
    J = get_jac(iter) #.cache.J
    u = iter.point.s ./ iter.point.y;
    start_advanced_timer(timer, "symmetric/form_system/update_indicies");
    #first_para_indicies, para_row_info = compute_indicies(J,u) # ONLY DO THIS ONCE???
    #kkt_solver.first_para_indicies = first_para_indicies
    first_para_indicies = kkt_solver.first_para_indicies
    para_row_info = kkt_solver.para_row_info
    update_indicies!(para_row_info, u)

    pause_advanced_timer(timer, "symmetric/form_system/update_indicies");

    J_new = J[first_para_indicies,:]
    #J_T_new = J_new'
    m_new = length(first_para_indicies)
    u_new = zeros(m_new) # the new diagonal elements
    for i = 1:m_new
        u_new[i] = para_row_info[i].u
    end

    if !all(u_new .>= 0)
        error("minimum(u_new) = $(minimum(u_new)) < 0.0")
    end

    start_advanced_timer(timer, "symmetric/form_system/M");
    U_new = spdiagm(u_new)
    n_var = size(J_new,2)
    M = [[get_lag_hess(iter) spzeros(n_var,m_new)]
        [J_new -U_new]];
    pause_advanced_timer(timer, "symmetric/form_system/M");

    pause_advanced_timer(timer, "symmetric/form_system");

    start_advanced_timer(timer, "symmetric/allocate");
    kkt_solver.Q = M;
    kkt_solver.factor_it = iter;
    kkt_solver.schur_diag = compute_schur_diag(iter)
    kkt_solver.true_x_diag = diag(M)[1:dim(iter)]
    kkt_solver.ready = :system_formed
    pause_advanced_timer(timer, "symmetric/allocate");
end

function factor_implementation!(kkt_solver::Clever_Symmetric_KKT_solver, timer::class_advanced_timer)
    num_var = dim(kkt_solver.factor_it)
    num_con_reduced = length(kkt_solver.para_row_info)
    @assert(size(kkt_solver.Q,2) == num_var + num_con_reduced)
    return ls_factor!(kkt_solver.ls_solver, kkt_solver.Q, num_var, num_con_reduced, timer)
end

function ls_solve(mat, solver, my_rhs::Array{Float64,1}, timer::class_advanced_timer, num_ref::Int64)
	# do iterative_refinement for num_ref iterations
	sol = zeros(length(my_rhs))
	for i = 1:num_ref
		if i == 1
			err = deepcopy(my_rhs)
		else
			err = my_rhs - vector_product(mat,sol)
		end
		sol += ls_solve(solver, err, timer)
		#@show i, norm(my_rhs - mat * sol)
	end
	return sol
end

function compute_direction_implementation!(kkt_solver::Clever_Symmetric_KKT_solver, timer::class_advanced_timer)
    start_advanced_timer(timer, "symmetric")
    y_org = get_y(kkt_solver.factor_it)
    s_org = get_s(kkt_solver.factor_it)
    rhs = kkt_solver.rhs;

    start_advanced_timer(timer, "symmetric/rhs_clever")
    symmetric_primal_rhs = rhs.primal_r + rhs.comp_r ./ y_org

    m = length(kkt_solver.para_row_info)
    clever_symmetric_primal_rhs = zeros(m)
    for i = 1:m
        for row in kkt_solver.para_row_info[i].ls
            j = row.ind
            clever_symmetric_primal_rhs[i] += row.g * symmetric_primal_rhs[j]
        end
    end

    clever_symmetric_rhs = [rhs.dual_r; clever_symmetric_primal_rhs];
    pause_advanced_timer(timer, "symmetric/rhs_clever")

    # DO ITERATIVE REFINEMENT!
    #@show kkt_solver.ready

    if false
        dir_x_and_y = ls_solve(kkt_solver.ls_solver, clever_symmetric_rhs, timer)
        #display(full(kkt_solver.Q))
        #@show norm(clever_symmetric_rhs - vector_product(kkt_solver.Q,dir_x_and_y))
        #@show norm(clever_symmetric_rhs - kkt_solver.Q * dir_x_and_y)
    else
        num_ref = 2
        dir_x_and_y = ls_solve(kkt_solver.Q, kkt_solver.ls_solver, clever_symmetric_rhs, timer, num_ref)
    end

    dir = kkt_solver.dir;
    dir.x = dir_x_and_y[1:length(rhs.dual_r)];
    v = dir_x_and_y[(1+length(rhs.dual_r)):end];

    start_advanced_timer(timer, "symmetric/dir_clever")
    y = zeros(length(dir.y))
    y_diag_vals = s_org ./ y_org;
    for j = 1:length(y)
        y[j] = y_diag_vals[j]^(-1.0) * symmetric_primal_rhs[j]
    end

    for i = 1:m
        row_group = kkt_solver.para_row_info[i]
        tmp = -(clever_symmetric_primal_rhs[i] + row_group.u * v[i])

        for row in row_group.ls
            j = row.ind
            y[j] += (row.u)^(-1.0) * row.ratio * tmp
        end
    end
    dir.y = y
    # dir.y = ...
    # dir.s = ( rhs.comp_r - dir.y .* s_org ) ./ y_org
    dir.s = eval_jac_prod(kkt_solver.factor_it, dir.x) - rhs.primal_r

    check_for_nan(dir)
    pause_advanced_timer(timer, "symmetric/dir_clever")

    start_advanced_timer(timer, "symmetric/kkt_err");
    update_kkt_error!(kkt_solver, Inf, timer)
    pause_advanced_timer(timer, "symmetric/kkt_err");
    pause_advanced_timer(timer, "symmetric")
end

function update_delta_vecs!(kkt_solver::Clever_Symmetric_KKT_solver, delta_x_vec::Array{Float64,1}, delta_s_vec::Array{Float64,1}, timer::class_advanced_timer)
    start_advanced_timer(timer, "symmetric")
    start_advanced_timer(timer, "symmetric/delta_vecs")
    kkt_solver.delta_x_vec = delta_x_vec
    kkt_solver.delta_s_vec = delta_s_vec

    if sum(abs.(delta_s_vec)) > 0.0
        error("not implemented")
    else
        if sum(abs.(delta_x_vec)) > 0.0
            n = length(delta_x_vec)
            m = size(kkt_solver.Q,2) - n
            kkt_solver.Q = kkt_solver.Q + spdiagm([rand(length(delta_x_vec)); zeros(m)])
            for i = 1:length(delta_x_vec)
                kkt_solver.Q[i,i] = kkt_solver.true_x_diag[i] + delta_x_vec[i]
            end
        end

    end

    kkt_solver.ready = :delta_updated
    pause_advanced_timer(timer, "symmetric/delta_vecs")
    pause_advanced_timer(timer, "symmetric")
end


text:
##### recombinations.jl #####

# In this file you will find all the functions regarding recombinations and crossovers.
# Recombinations are used specially for Evolution Strategies, while crossovers are used
# specially for Genetic Algorithms.

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

export crossover

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

# Recombinations
# ==============
function average(population ::Vector{T}) where {T <: Vector}
    obj = zeros(eltype(T), length(population[1]))
    l   = length(population)
    for i in 1:l
        obj += population[i]
    end
    return obj ./ l
end

function marriage(population ::Vector{T}) where {T <: Vector}
    s   = length(population)
    l   = length(population[1])
    obj = zeros(eltype(T), l)
    for i in 1:l
        obj[i] = population[rand(1:s)][i]
    end
    return obj
end

# Strategy recombinations
# =======================
function averageSigma1(ss ::Vector{S}) where {S <: Strategy}
    s = copy(ss[1])
    σ = 0.0
    l = length(ss)
    for i in 1:l
        σ += ss[i][:σ]
    end
    s[:σ] = σ / l
    return s
end

function averageSigmaN(ss ::Vector{S}) where {S <: Strategy}
    s = copy(ss[1])
    σ = zeros(length(ss[1][:σ]))
    l = length(ss)
    for i in 1:l
        σ += ss[i][:σ]
    end
    s[:σ] = σ ./ l
    return s
end

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

# ==================
# Genetic algorithms
# ==================

# Binary crossovers
# -----------------

# Single point crossover
function singlepoint(v1 ::T, v2 ::T) where {T <: AbstractVector}
    l   = length(v1)
    c1  = copy(v1)
    c2  = copy(v2)
    pos = rand(1:l)
    for i in pos:l
        vswap!(c1, c2, i)
    end
    return c1, c2
end

# Two point crossover
function twopoint(v1 ::T, v2 ::T) where {T <: AbstractVector}
    l  = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)
    from, to = rand(1:l, 2)
    from, to = from > to ? (to, from)  : (from, to)
    for i in from:to
        vswap!(c1, c2, i)
    end
    return c1, c2
end

# Uniform crossover
function uniform(v1 ::T, v2 ::T) where {T <: AbstractVector}
    l   = length(v1)
    c1  = copy(v1)
    c2  = copy(v2)
    for i in 1:l
        if rand(Bool)
            vswap!(c1, c2, i)
        end
    end
    return c1, c2
end

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

# Real valued crossovers
# ----------------------

# Discrete Crossover
function discrete(v1 ::T, v2 ::T) where {T <: AbstractVector}
    l    = length(v1)
    c1   = similar(v1)
    c2   = similar(v2)
    for i in 1:l
        @inbounds c1[i] = rand(Bool) ? v1[i] : v2[i]
        @inbounds c2[i] = rand(Bool) ? v2[i] : v1[i]
    end
    return c1, c2
end

# Weighted arithmetic mean
function waverage(w ::Vector{Float64})
    function wavexvr(v1 ::T, v2 ::T) where {T <: AbstractVector}
        c1 = (v1+v2) ./ w
        return c1, copy(c1)
    end
    return wavexvr
end

# Intermediate recombination
function intermediate(d ::Float64 = 0.0)
    function intermxvr(v1 ::T, v2 ::T) where {T <: AbstractVector}
        l  = length(v1)
        α  = (1.0+2d) * rand(l) .- d
        c1 = v2 .+ α .* (v1 - v2)
        α  = (1.0+2d) * rand(l) .- d
        c2 = v1 .+ α .* (v2 - v1)
        return c1, c2
    end
    return intermxvr
end

# Line recombination
function line(d ::Float64 = 0.0)
    function linexvr(v1 ::T, v2 ::T) where {T <: AbstractVector}
        α1, α2 = (1.0+2d) * rand(2) .- d
        c1 = v2 .+ α2 * (v1 - v2)
        c2 = v1 .+ α1 * (v2 - v1)
        return c1, c2
    end
    return linexvr
end

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

# Permutation crossovers
# ----------------------

# Partially mapped crossover
function pmx(v1 ::T, v2 ::T) where {T <: AbstractVector}
    s = length(v1)
    from, to = rand(1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = similar(v1)
    c2 = similar(v2)

    # Swap
    @inbounds c1[from:to] = v2[from:to]
    @inbounds c2[from:to] = v1[from:to]

    # Fill in from parents
    for i in vcat(1:from-1, to+1:s)
        # Check conflicting offspring
        in1 = inmap(v1[i], c1, from, to)
        if in1 == 0
            @inbounds c1[i] = v1[i]
        else
            tmpin = in1
            while tmpin > 0
                @inbounds tmpin = inmap(c2[in1], c1, from, to)
                in1   = tmpin > 0 ? tmpin : in1
            end
            @inbounds c1[i] = v1[in1]
        end

        in2 = inmap(v2[i], c2, from, to)
        if in2 == 0
            @inbounds c2[i] = v2[i]
        else
            tmpin = in2
            while tmpin > 0
                @inbounds tmpin = inmap(c1[in2], c2, from, to)
                in2   = tmpin > 0 ? tmpin : in2
            end
            @inbounds c2[i] = v2[in2]
        end
    end
    return c1, c2
end

# Order crossover
function ox1(v1 ::T, v2 ::T) where {T <: AbstractVector}
    s = length(v1)
    from, to = rand(1:s, 2)
    from, to = from > to ? (to, from)  : (from, to)
    c1 = zeros(v1)
    c2 = zeros(v2)
    # Swap
    @inbounds begin
        c1[from:to] = v2[from:to]
        c2[from:to] = v1[from:to]
    end
    # Fill in from parents
    k = to+1 > s ? 1 : to+1 #child1 index
    j = to+1 > s ? 1 : to+1 #child2 index
    for i in vcat(to+1:s,1:from-1)
        while in(v1[k],c1)
            k = k+1 > s ? 1 : k+1
        end
        @inbounds c1[i] = v1[k]
        while in(v2[j],c2)
            j = j+1 > s ? 1 : j+1
        end
        @inbounds c2[i] = v2[j]
    end
    return c1, c2
end

# Cycle crossover
function cx(v1 ::T, v2 ::T) where {T <: AbstractVector}
    s  = length(v1)
    c1 = zeros(v1)
    c2 = zeros(v2)

    f1 = true #switch
    k  = 1
    while k > 0
        idx = k
        if f1
            #cycle from v1
            while c1[idx] == zero(T)
                @inbounds c1[idx] = v1[idx]
                @inbounds c2[idx] = v2[idx]
                idx = inmap(v2[idx],v1,1,s)
            end
        else
            #cycle from v2
            while c2[idx] == zero(T)
                @inbounds c1[idx] = v2[idx]
                @inbounds c2[idx] = v1[idx]
                idx = inmap(v1[idx],v2,1,s)
            end
        end
        f1 $= true
        k = inmap(zero(T),c2,1,s)
    end
    return c1,c2
end

# Order-based crossover
function ox2(v1 ::T, v2 ::T) where {T <: AbstractVector}
    s  = length(v1)
    c1 = copy(v1)
    c2 = copy(v2)

    for i in 1:s
        if rand(Bool)
            @inbounds begin
                idx1 = inmap(v2[i],v1,1,s)
                idx2 = inmap(v1[i],v2,1,s)
                c1[idx1] = zero(T)
                c2[idx2] = zero(T)
            end
        end
    end

    for i in 1:s
        if !in(v2[i],c1)
            tmpin     = inmap(zero(T),c1,1,s)
            @inbounds c1[tmpin] = v2[i]
        end
        if !in(v1[i],c2)
            tmpin     = inmap(zero(T),c2,1,s)
            @inbounds c2[tmpin] = v1[i]
        end
    end
    return c1,c2
end

# Position-based crossover
function pos(v1 ::T, v2 ::T) where {T <: AbstractVector}
    s  = length(v1)
    c1 = zeros(v1)
    c2 = zeros(v2)

    for i in 1:s
        if rand(Bool)
            @inbounds begin
                c1[i] = v2[i]
                c2[i] = v1[i]
            end
        end
    end

    for i in 1:s
        if !in(v1[i],c1)
            tmpin     = inmap(zero(T),c1,1,s)
            @inbounds c1[tmpin] = v1[i]
        end
        if !in(v2[i],c2)
            tmpin     = inmap(zero(T),c2,1,s)
            @inbounds c2[tmpin] = v2[i]
        end
    end
    return c1, c2
end

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

"""
    crossover(gene1 ::T, gene2 ::T) where {T <: IntegerGene}

Crosses two `IntegerGene` genes according to the crossover function chosen.
"""
function crossover(gene1 ::T, gene2 ::T) where {T <: IntegerGene}
    return crossover(gene1.value, gene2.value)
end

"""
    crossover(gene1 ::T, gene2 ::T) where {T <: FloatGene}

Crosses two `FloatGene` genes according to the crossover function chosen.
"""
function crossover(gene1 ::T, gene2 ::T) where {T <: FloatGene}
    return crossover(gene1.value, gene2.value)
end

"""
    crossover(chromo1 ::T, chromo2 ::T) where {T <: Vector{<:AbstractGene}}

`chromo1` and `chromo2` are two vectors of genes. Crosses each entry of both chromossomes according to the crossover function chosen.
"""
function crossover(chromo1 ::T, chromo2 ::T) where {T <: Individual}
    c1 = deepcopy(chromo1)
    c2 = deepcopy(chromo2)
    for i in 1:length(chromo1)
        @inbounds c1[i].value, c2[i].value = crossover(chromo1[i], chromo2[i])
    end
    return c1, c2
end

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

# Utils
# =====
function vswap!(v1 ::T, v2 ::T, idx ::Int) where {T <: AbstractVector}
    @inbounds begin
        val     = v1[idx]
        v1[idx] = v2[idx]
        v2[idx] = val
    end
end

function inmap(v ::T, c ::Vector{T}, from ::Int, to ::Int) where {T}
    exists = 0
    for j in from:to
        if exists == 0 && v == c[j]
            exists = j
        end
    end
    return exists
end


text:
__precompile__()

module Pincus

using Base.Cartesian, FastGaussQuadrature, Sobol#, CLArrays, GPUArrays

export  gauss_hermite_tup,
        uniboxmaximize,
        threaded_uniboxmaximize,
        maximize,
        threaded_maximize,
        sobol_vec

include("opt.jl")
include("sobol.jl")

end # module


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using JSON
export LocalSearchAlgorithm, NeighborhoodAlgorithm, SearchGraph, find_neighborhood, push_neighborhood!

abstract type LocalSearchAlgorithm end
abstract type NeighborhoodAlgorithm end

mutable struct SearchGraph{T} <: Index
    db::Vector{T}
    recall::Float64
    k::Int
    links::Vector{Vector{Int32}}
    search_algo::LocalSearchAlgorithm
    neighborhood_algo::NeighborhoodAlgorithm
    verbose::Bool
end

@enum VertexSearchState begin
    UNKNOWN = 0
    VISITED = 1
    EXPLORED = 2
end

function fit(::Type{SearchGraph}, dist::Function, dataset::AbstractVector{T}; recall=0.9, k=10, search_algo=BeamSearch(), neighborhood_algo=LogSatNeighborhood(1.1), verbose=true) where T
    links = Vector{Int32}[]
    index = SearchGraph(T[], recall, k, links, search_algo, neighborhood_algo, verbose)
    for item in dataset
        push!(index, dist, item)
    end

    index
end

include("opt.jl")

## neighborhoods
include("neighborhood/fixedneighborhood.jl")
include("neighborhood/logneighborhood.jl")
include("neighborhood/logsatneighborhood.jl")
include("neighborhood/gallopingneighborhood.jl")
include("neighborhood/essencialneighborhood.jl")
include("neighborhood/satneighborhood.jl")
include("neighborhood/galsatneighborhood.jl")
include("neighborhood/vorneighborhood.jl")

## search algorithms
include("ihc.jl")
include("tihc.jl")
include("beamsearch.jl")

### Basic operations on the index
const OPTIMIZE_LOGBASE = 10
const OPTIMIZE_LOGBASE_STARTING = 4

"""
    find_neighborhood(index::SearchGraph{T}, dist::Function, item)

Searches for `item` neighborhood in the index, i.e., if `item` were in the index whose items should be
its neighbors (intenal function)
"""
function find_neighborhood(index::SearchGraph, dist::Function, item)
    n = length(index.db)
    n == 0 && return (KnnResult(1), Int32[])
    neighborhood(index.neighborhood_algo, index, dist, item)
end

"""
    push_neighborhood!(index::SearchGraph, item, L::AbstractVector{Int32}, n::Int)

Inserts the object `item` into the index, i.e., creates an edge from items listed in L and the
vertex created for ìtem` (internal function)
"""
function push_neighborhood!(index::SearchGraph{T}, item::T, L::AbstractVector{Int32}, n::Int) where T
    for objID in L
        push!(index.links[objID], 1+n)
    end

    push!(index.links, L)
    push!(index.db, item)
end

"""
    push!(index::SearchGraph, dist::Function, item)

Inserts `item` into the index.
"""
function push!(index::SearchGraph, dist::Function, item) where T
    knn, neighbors = find_neighborhood(index, dist, item)
    n = length(index.db)
    push_neighborhood!(index, item, neighbors, n)
    n += 1
    if n > OPTIMIZE_LOGBASE_STARTING
        k = ceil(Int, log(OPTIMIZE_LOGBASE, 1+n))
        k1 = ceil(Int, log(OPTIMIZE_LOGBASE, 2+n))
        k != k1 && optimize!(index, dist, recall=index.recall)
    end

    if index.verbose && length(index.db) % 5000 == 0
        println(stderr, "added n=$(length(index.db)), neighborhood=$(length(neighbors)), $(now())")
    end
    knn, neighbors
end

const EMPTY_INT_VECTOR = Int[]

"""
    search(index::SearchGraph, dist::Function, q, res::KnnResult; hints)

Solves the specified query `res` for the query object `q`.
If hints is given then these vertices will be used as starting poiints for the search process.
"""
function search(index::SearchGraph, dist::Function, q, res::KnnResult; hints=EMPTY_INT_VECTOR)
    length(index.db) == 0 && return res
    navigation_state = Dict{Int,VertexSearchState}()
    sizehint!(navigation_state, maxlength(res))
    search(index.search_algo, index, dist, q, res, navigation_state, hints)
end

"""
    optimize!(index::SearchGraph, dist::Function;
        recall=0.9,
        k=10,
        num_queries=128,
        perf=nothing,
        tol::Float64=0.01,
        maxiters::Int=3,
        probes::Int=0)

Optimizes the index for the specified kind of queries.
"""
function optimize!(index::SearchGraph{T}, dist::Function;
    recall=0.9,
    k=10,
    num_queries=128,
    perf=nothing,
    tol::Float64=0.01,
    maxiters::Int=3,
    probes::Int=0) where T
    if perf === nothing
        perf = Performance(index.db, dist; expected_k=k, num_queries=num_queries)
    end
    optimize!(index.search_algo, index, dist, recall, perf, tol=tol, maxiters=3, probes=probes)
end


text:
module PPM

function write(pix)
  w, h = size(pix)

  println("P3")
  print(w)
  print(" ")
  println(h)
  println(0xffff)

  for iw = 1:w
    for ih = 1:h
      print(pix[iw,ih][1], " ", pix[iw,ih][2], " ", pix[iw,ih][3])
      print("  ")
    end
    println()
  end
end

end

max_stars_repo_name: dominichamon/juliajulia

max_issues_repo_name: dominichamon/juliajulia

max_forks_repo_name: dominichamon/juliajulia


text:
using BandedMatrices, Test

@testset "QR tests" begin
    for T in (Float64,ComplexF64,Float32,ComplexF32)
        A=brand(T,10,10,3,2)
        Q,R=qr(A)
        @test Matrix(Q)*Matrix(R) ≈ A
        b=rand(T,10)
        @test mul!(similar(b),Q,mul!(similar(b),Q',b)) ≈ b
        for j=1:size(A,2)
            @test Q' * A[:,j] ≈ R[:,j]
        end
        A=brand(T,14,10,3,2)

        Q,R=qr(A)

        for k=1:size(Q,1),j=1:size(Q,2)
            @test Q[k,j] ≈ Matrix(Q)[k,j]
        end

        @test Matrix(Q)*Matrix(R) ≈ A

        A=brand(T,10,14,3,2)

        Q,R=qr(A)

        for k=1:size(Q,1),j=1:size(Q,2)
            @test Q[k,j] ≈ Matrix(Q)[k,j]
        end

        @test Matrix(Q)*Matrix(R) ≈ A
        A=brand(T,100,100,3,4)
        Q,R=qr(A)
        b=rand(T,100)
        @test R\(Q'*b) ≈ Matrix(A)\b
    end

    A=brand(10,10,3,2)
    b=rand(ComplexF64,10)
    Q,R=qr(A)
    @test R\(Q'*b) ≈ Matrix(A)\b

    A=brand(ComplexF64,10,10,3,2)
    b=rand(10)
    Q,R=qr(A)
    @test R\(Q'*b) ≈ Matrix(A)\b
end