text:
include("meta.jl")

#@enum Result res_mem res_madd res_madd_clamp res_zero res_one res_m_one res_two res_1_fifth res_m_6_fifths res_exp2
@constants result_selections begin
	res_maxval # 0
	res_madd   # 1
	res_mem    # 2
	res_ext    # 3
	res_one    # even
	res_m_one  # odd
	res_two    # even
	res_m_two  # odd
	res_lfsr   # even
	res_eps    # odd
	res_1_fifth
	res_m_6_fifths
	res_zero
	res_exp2
end

const res_const = res_maxval

const const_sel_bits = 5
const result_sel_bits = 4
@flags control_flags 4 begin
# Data path control flags
	# register write enable bits
	T_WE
	F1_WE
	F2_WE
	SIGN_WE   # store sign(f2) into dp.positive
	OVERFLOW_WE # store overflow flag from madd result
	OUTPUT_WE # should probably not be handled by data path
	BIAS_WE
	# other control flags
	C_CLAMP
	C_STORED_SIGN
	C_TRUNC
	C_DITHER

# Non-data path control flags
	C_M_WE
	C_CONDITIONAL_M_WE
	C_MEM_ADDR
	C_INST_DONE
end

const res_madd_clamp = res_madd | C_CLAMP;
const res_stored_sign = res_one | C_STORED_SIGN;

const A_WE = T_WE | F1_WE | F2_WE | C_TRUNC

const MEM_READ = C_MEM_ADDR
const IMM_WE = C_M_WE
const MEM_WE = C_M_WE | C_MEM_ADDR

struct DPControl
	reg_we::Int
	res::Int
end
DPControl(reg_we, res) = DPControl(reg_we, res, 0)

const IC_N_ZERO = -13  # When n == IC_N_ZERO in IContext, the result is scale = 0

struct IContext
	addr::Int
	n::Int
end
IContext(addr) = IContext(addr, 0)
IContext() = IContext(0)

struct CycleControl
	cdp::DPControl
	cn::Int
end

control(flags::Int, res::Int) = CycleControl(DPControl(flags, res), flags)
get_res(c::Int) = c & ((1 << result_sel_bits) - 1)
get_flags(c::Int) = c & (-1 << result_sel_bits)

get_bits(c::CycleControl) = c.cdp.res | c.cn
get_res(c::CycleControl) = get_res(c.cdp.res)
get_flags(c::CycleControl) = get_flags(get_bits(c))

@enum Opcode begin
	inst_nop

	inst_acc_imm_to_mem
	inst_macc_to_mem
	inst_wadd
	inst_lpfilter

	inst_pulse_exp2
	inst_pulse
	inst_triangle
	inst_sina2
	inst_sin

	inst_lfsr

	inst_output
	inst_output_a
	inst_dither_update

	inst_imm
	inst_acc_to_imm
end

struct Instruction
	op::Opcode
	context::IContext
end
Instruction(op::Opcode) = Instruction(op, IContext())

const microcode = [
	inst_nop => [control(0, res_mem)],
	inst_imm => [control(A_WE, res_mem)], # a = imm
	inst_acc_to_imm => [ # a = (imm += 2^-n*a + dither)
		control(F1_WE, res_exp2),
		control(T_WE,  res_mem),
		control(A_WE | IMM_WE | OVERFLOW_WE | C_DITHER, res_madd)
	],
	inst_acc_imm_to_mem => [ # a = (mem += 2^-n*imm + dither)
		control(F1_WE, res_exp2),
		control(F2_WE, res_mem),
		control(T_WE | MEM_READ, res_mem),
		control(A_WE | MEM_WE | OVERFLOW_WE | C_DITHER, res_madd)
	],
	inst_wadd => [ # a = 4^-n*(imm*a) + mem (n can be < 0)
		# f2 = imm*a
		control(T_WE, res_zero),
		control(F1_WE, res_mem),
		control(F2_WE, res_madd),
		# f2 = 2^-n*f2
		control(F1_WE, res_exp2),
		control(F2_WE, res_madd),
		# a = 2^-n*f2 + mem
		control(T_WE | MEM_READ, res_mem),
		control(A_WE, res_madd)
	],
	inst_macc_to_mem => [ # mem = clamp(imm*a + mem) -- keeps value of a!
		control(F1_WE, res_mem),
		control(T_WE | MEM_READ, res_mem),
		control(MEM_WE, res_madd_clamp),
		control(F1_WE, res_one), # restore a from f2...
		control(T_WE, res_zero), #
		control(A_WE, res_madd)  #
	],
	inst_lpfilter =>  [ # a = mem = clamp(mem + 2^-n*(imm - mem) + dither)
		# f2 = imm-mem
		control(T_WE, res_mem),
		control(F1_WE, res_m_one),
		control(F2_WE | MEM_READ, res_mem),
		control(F2_WE, res_madd_clamp),
		# mem += 2^-n*f2 + dither (clamp)
		control(T_WE | MEM_READ, res_mem),
		control(F1_WE, res_exp2),
		control(A_WE | MEM_WE, res_madd_clamp | C_DITHER)
	],
	inst_sina2 => [ # a = -sina2(a)
		# a = 2a + 1/2
		control(F1_WE | SIGN_WE, res_two),
		control(T_WE, res_one),
		control(F1_WE | F2_WE | C_TRUNC, res_madd),
		# a = a*a - 1
		control(T_WE, res_m_one),
		control(F2_WE, res_madd),
		# a *= sign(a_in)
		control(F1_WE, res_stored_sign),
		control(T_WE, res_zero),
		control(A_WE, res_madd_clamp),
	],
	inst_sin => [ # a = sina4(a)
		# a = 2a + 1/2
		control(F1_WE | SIGN_WE, res_two),
		control(T_WE, res_one),
		control(F1_WE | F2_WE | C_TRUNC, res_madd),
		# f1 = a*a
		control(T_WE, res_zero),
		control(F1_WE, res_madd),
		# f2 = -6/5 + 1/5*f1
		control(T_WE, res_m_6_fifths),
		control(F2_WE, res_1_fifth),
		control(F2_WE, res_madd),
		# a = M + f1*f2
		control(T_WE, res_one),
		control(F2_WE, res_madd),
		# a *= sign(a_in)
		control(F1_WE, res_stored_sign),
		control(T_WE, res_zero),
		control(A_WE, res_madd_clamp),
	],
	inst_output => [ # output imm; a = imm = 0
		control(OUTPUT_WE, res_mem),
		control(A_WE | IMM_WE, res_zero)
	],
	inst_dither_update => [ # bias = bit_reverse(imm); a = (imm += epsilon)
		control(T_WE | BIAS_WE, res_mem),
		control(F1_WE, res_one),
		control(F2_WE, res_eps),
		control(A_WE | IMM_WE, res_madd)
	],
	inst_output_a => [ # output a
		control(F1_WE | F2_WE, res_zero),
		control(A_WE | OUTPUT_WE, res_madd)
	],
	inst_pulse => [ # a = a >= imm ? maxval : minval
		control(F1_WE, res_m_one),
		control(F2_WE, res_mem),
		control(F2_WE, res_madd), # f2 = a - imm
		control(SIGN_WE, res_zero), # s = sign(a - imm), TODO: Merge with previous cycle when sign_we reads from result!
		control(A_WE, res_maxval | C_STORED_SIGN)
	],
	inst_pulse_exp2 => [ # a = a >= -1 + 2^-n ? maxval : minval
		# Prepare
		control(F1_WE | F2_WE, res_one),
		control(T_WE, res_madd),
		# Same as inst_pulse except f2 = 2^-n instead of f2 = imm
		# TODO: Merge with inst_pulse when 2^-n can represent imm
		control(F1_WE, res_m_one),
		control(F2_WE, res_exp2),
		control(F2_WE, res_madd), # f2 = a - imm
		control(SIGN_WE, res_zero), # s = sign(a - imm), TODO: Merge with previous cycle when sign_we reads from result!
		control(A_WE, res_maxval | C_STORED_SIGN)
	],
	inst_triangle => [ # a = triangle(a)
		# a += 2^-n -- use to get a += 0.5
		control(F1_WE, res_one),
		control(F2_WE, res_exp2),
		control(A_WE, res_madd | C_TRUNC),

		control(T_WE | SIGN_WE, res_m_one),
		control(F1_WE, res_two | C_STORED_SIGN),
		control(A_WE, res_madd_clamp)
	],
	inst_lfsr => [ # a = imm = lfsr(a)
		control(F1_WE, res_two),
		control(F2_WE, res_mem),
		control(T_WE, res_lfsr),
		control(A_WE | IMM_WE | C_CONDITIONAL_M_WE, res_madd)
	],
]
const opcode_programs = Dict(microcode)

const instructions = [
	"nop"            => Instruction(inst_nop,            IContext(0,  IC_N_ZERO)),

	"enable"         => Instruction(inst_nop,            IContext(0,  IC_N_ZERO+2)),
	"disable"        => Instruction(inst_nop,            IContext(0,  IC_N_ZERO+3)),

	"phase_update"   => Instruction(inst_acc_imm_to_mem, IContext(0,  IC_N_ZERO)),
	"contribute"     => Instruction(inst_macc_to_mem,    IContext(0,  IC_N_ZERO)),
	"madd_scale2"    => Instruction(inst_wadd,           IContext(0,  IC_N_ZERO)),
	"approach"       => Instruction(inst_lpfilter,       IContext(0,  IC_N_ZERO)),

	"square"         => Instruction(inst_pulse_exp2,     IContext(0,  0)),
	"pulse"          => Instruction(inst_pulse_exp2,     IContext(0,  IC_N_ZERO)),
	"pulse_imm"      => Instruction(inst_pulse,          IContext(0,  IC_N_ZERO)),
	"sawtooth"       => Instruction(inst_nop,            IContext(0,  IC_N_ZERO)), # synonym for nop
	"triangle"       => Instruction(inst_triangle,       IContext(0, -1)),
	"sina2"          => Instruction(inst_sina2,          IContext(0,  IC_N_ZERO)),
	"sin"            => Instruction(inst_sin,            IContext(0,  IC_N_ZERO)),

	"noise_update"   => Instruction(inst_lfsr,           IContext(0,  IC_N_ZERO)),

	"output"         => Instruction(inst_output,         IContext(0,  IC_N_ZERO)),
	"output_left"    => Instruction(inst_output,         IContext(0,  IC_N_ZERO)),
	"output_right"   => Instruction(inst_output,         IContext(0,  IC_N_ZERO+1)),
	"output_a"       => Instruction(inst_output_a,       IContext(0,  IC_N_ZERO)),
	"output_a_left"  => Instruction(inst_output_a,       IContext(0,  IC_N_ZERO)),
	"output_a_right" => Instruction(inst_output_a,       IContext(0,  IC_N_ZERO+1)),

	"loop_update"    => Instruction(inst_dither_update,  IContext(0,  IC_N_ZERO)),
]

max_stars_repo_name: toivoh/fmicrosynth

max_issues_repo_name: toivoh/fmicrosynth

max_forks_repo_name: toivoh/fmicrosynth


text:
@testset "System Tests" begin
    @testset "Composition System" begin
        @var x y a b

        f = System([y^2 + 2x + 3, x - 1])
        g = System([x + y * a, x - 1 * b], parameters = [a, b])
        @test parameters(g ∘ f) == [a, b]
        @test nparameters(g ∘ f) == 2

        @test parameters(g ∘ g) == [a, b]
        @test nparameters(g ∘ g) == 2

        @test parameters(f ∘ g) == [a, b]
        @test nparameters(f ∘ g) == 2

        @test parameters(f ∘ f) == []
        @test nparameters(f ∘ f) == 0

        @test parameters(f ∘ g ∘ f) == [a, b]
        @test nparameters(f ∘ g ∘ f) == 2
    end
end

ext: jl


text:
#TODO: find julia native alternative for it
#=
    This script just run all the other scripts in the correct order.
=#

## ------------------------------------------------------------------------
using ArgParse

# function check_clear_args(args)
#     for arg in split(args, ",")
#         !(arg in ["all", "base", "maxent_ep", "fva_pp", "cache"]) && return false
#     end
#     return true
# end

set = ArgParseSettings()
@add_arg_table! set begin
    "--force-pull"
        help = "make a sync with the remote"
        action = :store_true
    "--install"
        help = "run an installation script"
        action = :store_true
    "--clear-cache"
        help = "clear the cache folder"
        action = :store_true
    "--clear-fva-models"
        help = "clear the fva_pp_models folder"
        action = :store_true

    # "--dry-run"
    #     help = "run without consecuences, just printing"
    #     action = :store_true
    # "--run", "-r", "-x"
    #     help = "possible values: \"all\" (run all the scripts),  " *
    #                              "\"base\" (run only the base model scripts), " *
    #                              "\"none\" (run nothing)"
    #     default = "base"
    #     required = false
    #     range_tester = (x -> x in ["all", "base", "none"])
    # "--clear", "-c"
    #     help = "possible values: \"raw\" (clear raw data folder), " *
    #                             "\"all\" (clear all the scripts targets), " *
    #                             "\"base\" (clear only the base model scripts targets), " *
    #                             "\"maxent_ep\" (clear only the maxent_ep bundles), " *
    #                             "\"fva_pp\" (clear only the fva preprocess models), " *
    #                             "\"cache\" (clear the cache forder). " *
    #                             "You can pass several using comma Ex: --clear=cache,maxent"
    #     required = false
    #     range_tester = check_clear_args
end

if isinteractive()
    # Dev vals
    force_pull_flag = false
    install_flag = false
    clear_cache_flag = false
    clear_fva_models_flag = false
else
    parsed_args = parse_args(set)
    force_pull_flag = parsed_args["force-pull"]
    install_flag = parsed_args["install"]
    clear_cache_flag = parsed_args["clear-cache"]
    clear_fva_models_flag = parsed_args["clear-fva-models"]
    # clear_args = parsed_args["clear"]
    # clear_args = isnothing(clear_args) ? nothing : split(clear_args, ",")
end

## ------------------------------------------------------------------------
using Pkg
try import DrWatson
catch
    import Pkg
    pkg"add DrWatson"
end
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Human1")

## ------------------------------------------------------------------------
function force_pull()
    here = pwd()
    cd(@__DIR__)
    run(`git reset -q HEAD -- .`)
    run(` git checkout -- .`)
    try;
        run(`git pull`)
    catch err
        @warn string("fail to pull: error ", err)
    end
    cd(here)
end
force_pull_flag && force_pull()

## ------------------------------------------------------------------------
if install_flag
    # sync
    force_pull()
    # Install unregistered packages
    try
        pkg"rm Chemostat"
        pkg"rm UtilsJL"
    catch; end
    pkg"add https://github.com/josePereiro/UtilsJL.git#master"
    pkg"add https://github.com/josePereiro/Chemostat#adbeb2f"
    pkg"instantiate"
    pkg"build"
    pkg"test Chemostat"
end

## ------------------------------------------------------------------------
import Chemostat_Human1
const H1 = Chemostat_Human1
cd(H1.PROJ_ROOT)

## ------------------------------------------------------------------------
# clear stuff
if clear_cache_flag
    if isdir(H1.CACHE_DIR)
        rm(H1.CACHE_DIR, force = true, recursive = true)
        @info string(relpath(H1.CACHE_DIR), " deleted!!!")
    end
end

if clear_fva_models_flag
    if isdir(H1.FVA_PP_MODELS_DIR)
        rm(H1.FVA_PP_MODELS_DIR, force = true, recursive = true)
        @info string(relpath(H1.FVA_PP_MODELS_DIR), " deleted!!!")
    end
end


text:
using CalibrationErrors

@testset "Quadratic: Two-dimensional example" begin
    skce = QuadraticUnbiasedSKCE(UniformScalingKernel(4, SquaredExponentialKernel(2)))

    # only two predictions, i.e., one term in the estimator
    @test calibrationerror(skce, ([1 0; 0 1], [1, 2])) ≈ 0
    @test calibrationerror(skce, ([1 0; 0 1], [1, 1])) ≈ 0
    @test calibrationerror(skce, ([1 0; 0 1], [2, 1])) ≈ -8 * exp(-4)
    @test calibrationerror(skce, ([1 0; 0 1], [2, 2])) ≈ 0
end

@testset "Linear: Two-dimensional example" begin
    skce = LinearUnbiasedSKCE(UniformScalingKernel(4, SquaredExponentialKernel(2)))

    # only two predictions, i.e., one term in the estimator
    @test calibrationerror(skce, ([1 0; 0 1], [1, 2])) ≈ 0
    @test calibrationerror(skce, ([1 0; 0 1], [1, 1])) ≈ 0
    @test calibrationerror(skce, ([1 0; 0 1], [2, 1])) ≈ -8 * exp(-4)
    @test calibrationerror(skce, ([1 0; 0 1], [2, 2])) ≈ 0
end


text:
using HDF5, JLD
using PyCall
@pyimport matplotlib2tikz

include("CenterOfMass.jl")
include("setup_parameters.jl")

save_plots = false
do_belief = false

plt[:close]("all")
pygui(!save_plots)

folder = "multi_robot"

type TrialData
  belief::Histogram
  estimate
  error
  applied_f
  f_hat
  csqmi
  robots_measurement
  robots
end

@load "$(folder)/data.25_iter.0_4" data_array csqmi_ratios thetas actuator_limit

forces = map(x->x.applied_f, data_array)
num_measurement = map(x->length(x.robots_measurement), data_array)
num_robot = map(x->length(x.robots), data_array)
actuator_limits = actuator_limit * num_measurement
limits = forces - actuator_limits .>-1e-6
partial_limits = forces - actuator_limit .>-1e-6

#######################
# plot normalized error
#######################
println("normalized error")

figure()
errors = map(x->x.error, data_array)
mean_errors = mean(errors, 2)

colors = ["b", "g", "r", "c"]

indices = collect(1:size(mean_errors,3))
for ii = 1:size(mean_errors, 1)
#for ii = 1:3
  plot(indices, mean_errors[ii, 1, :][:], color = colors[ii], linewidth = 2.0)
end
legend(map(string, csqmi_ratios))

cutoff = 0.15
for ii = 1:size(errors, 1)
  n = 2
  stds = std(errors[ii,:,:], 2)[:]/sqrt(size(errors,2))
  es = mean_errors[ii,1,:][:]
  fill([indices;reverse(indices)], [es+stds;reverse(es-stds)], color = colors[ii],
    alpha=0.2, linewidth=0.0)
end
#plot([0;25],[cutoff;cutoff],linestyle="--", color="k")
#plot([9;9],[0.0;cutoff],linestyle="--", color="k")
#plot([13;13],[0.0;cutoff],linestyle="--", color="k")
#plot([17;17],[0.0;cutoff],linestyle="--", color="k")

if false
for jj = 1:size(errors, 1)
  for ii = 1:size(errors, 2)
    plot(1:size(errors,3), errors[jj, ii,:][:], color = colors[jj], alpha = 0.2,
    linewidth = 1.0)

    #robot_added_inds = find(diff(num_robot[jj, ii, :][:]))+1
    #scatter(robot_added_inds, errors[jj, ii, robot_added_inds], color =
    #colors[jj], alpha = 0.2, marker = "1")

    #for kk = 2:size(errors, 3)
      #if
        #scatter(kk, errors[jj, ii, kk], color = colors[jj], alpha = 0.2)
      #end
    #end
  end
end
end

#ylim(0.0, 0.8)
xlim(0.0, size(errors, 3))

xlabel("Iteration")
ylabel("Normalized error")

matplotlib2tikz.save("$(folder)/convergence.tex", figureheight="\\figureheight",
  figurewidth="\\figurewidth")

#######################
# plot number of robots
#######################
println("number robots")

figure()

mean_robots = mean(num_robot, 2)
for ii = 1:size(mean_robots, 1)
  plot(indices, mean_robots[ii, 1, :][:], color = colors[ii], linewidth = 2.0)
end

legend(map(string, csqmi_ratios))

for ii = 1:size(num_robot, 2)
  for jj = 1:size(num_robot, 1)
    plot(1:size(errors,3), num_robot[jj, ii,:][:], color = colors[jj], alpha = 0.2,
    linewidth = 1.0)
  end
end

xlabel("Iteration")
ylabel("Num. robots")

matplotlib2tikz.save("$(folder)/num_robots.tex", figureheight="\\figureheight",
  figurewidth="\\figurewidth")

#############
# plot belief
#############
println("belief")

if save_plots && do_belief
  for (ii, csqmi_ratio) = enumerate(csqmi_ratios)
    for jj = 1:size(errors,2)
      trial_folder = "$(folder)/$(csqmi_ratio)_$(jj)"
      mkpath(trial_folder)

      for kk = 1:size(errors,3)
        plt[:close]("all")

        println("belief: ratio=$(csqmi_ratio), trial=$(jj), iter=$(kk)")
        data = data_array[ii,jj,kk]
        robots = data.robots
        belief = data.belief
        robots_measurement = data.robots_measurement
        theta = thetas[:,jj]
        com_p = theta[1:2]
        mass = theta[3]

        figure()
        plot_attachment_points(attachment_ps)
        plot_occupied_points(attachment_ps[setdiff(robots, robots_measurement)])
        cloud = to_cloud(belief, interior_q)
        scaled_ps = (15*cloud[4,:]'/maximum(cloud[4,:]))
        scatter3D(cloud[1,:]', cloud[2,:]', cloud[3,:]', "z", scaled_ps.^2, "purple", alpha = 0.5)

        scatter3D([com_p[1]], [com_p[2]], [mass], "z", 200, "red", marker="*", alpha=1)

        plot_measurement_points(attachment_ps[robots_measurement])

        if kk > 1
          old_robots = data_array[ii,jj,kk-1].robots
          if length(robots) > length(old_robots)
            plot_new_point(attachment_ps[setdiff(robots, old_robots)[1]])
          end
        end
        fill3d(hcat(circle_ps...), alpha=0.2)

        xlabel("X")
        ylabel("Y")
        zlabel("M")
        #fig[:axes][1][:get_yaxis]()[:set_visible](false)
        #fig[:axes][1][:get_xaxis]()[:set_visible](false)
        #axis("off")
        xlim(-1.0,1.0)
        ylim(-1.0,1.0)
        zlim(0.0,1.405)
        savefig("$(trial_folder)/belief_$(kk).png", pad_inches=0.01, bbox_inches="tight")
      end
    end
  end
end

pygui(true)

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
velocity(r::AbstractArray, t::Int, δt::Int) = (r[t+δt, :] .- r[t, :]) ./ δt
c(v1::AbstractVector, v2::AbstractVector) = v1 ⋅ v2

function vacf(
    df::DataFrame,
    id::Symbol,
    x::Symbol,
    y::Symbol,
    max_τ::Int,
    δt::Int,
)
    result = DataFrame(n = Int[], τ = Int[], value = Float64[])
    @inbounds for n = 1:maximum(df[!, id])
        data = extract(df, n, id, [x, y])
        T = size(data, 1)
        T < max_τ ? maxt = T : maxt = max_τ
        @inbounds for τ = 0:maxt-δt-1
            s = 0.0
            @simd for t = 1:(T-δt-τ)
                vₜ = velocity(data, t, δt)
                v2 = velocity(data, t + τ, δt)
                s += c(vₜ, v2)
            end
            s /= (T - δt - τ)
            push!(result, [n, τ, s])
        end
    end
    result
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
using Test
using StatsModelComparisons
using StanSample, StatsFuns
using Printf
using JSON

@testset "Arsenic"  begin

    ProjDir = @__DIR__

    #=
    if haskey(ENV, "JULIA_CMDSTAN_HOME")
        include(joinpath(ProjDir, "test_demo_wells.jl"))
    else
      println("\nJULIA_CMDSTAN_HOME not set. Skipping tests")
    end
    =#

    include(joinpath(ProjDir, "cvit.jl"))

    # Data
    data = JSON.parsefile(joinpath(ProjDir, "wells.data.json"))
    y = Float64.(data["switched"])
    x = Float64[data["arsenic"]  data["dist"]]
    n, m = size(x)

    # Model
    model_str = read(open(joinpath(ProjDir, "arsenic_logistic.stan")), String)
    sm1 = SampleModel("arsenic_logistic", model_str)

    data1 = (p = m, N = n, y = Int.(y), x = x)
    # Fit the model in Stan
    rc1 = stan_sample(sm1; data=data1)
    if success(rc1)
        nt1 = read_samples(sm1)

        # Compute LOO and standard error
        log_lik = nt1.log_lik'
        loo, loos, pk = psisloo(log_lik)
        elpd_loo = sum(loos)
        se_elpd_loo = std(loos) * sqrt(n)

        @test elpd_loo ≈ -1968.3 atol=2.0
        @test se_elpd_loo ≈ 15.5 atol=0.5
        @test all(pk .< 0.5)
    end
    println()

    # Fit a second model, using log(arsenic) instead of arsenic
    x2 = Float64[log.(data["arsenic"])  data["dist"]]

    # Model
    data2 = (p = m, N = n, y = Int.(y), x = x2)
    # Fit the model in Stan
    rc2 = stan_sample(sm1; data=data2)

    if success(rc2)
        nt2 = read_samples(sm1)
        # Compute LOO and standard error
        log_lik = nt2.log_lik'
        loo2, loos2, pk2 = psisloo(log_lik)
        elpd_loo = sum(loos2)
        se_elpd_loo = std(loos2) * sqrt(n)

        @test elpd_loo ≈ -1952.1 atol=2.0
        @test se_elpd_loo ≈ 16.2 atol=0.5
        @test all(pk .< 0.5)
    end

    if success(rc1) && success(rc2)
        ## Compare the models
        loodiff = loos - loos2

        @test sum(loodiff) ≈ -16.3 atol=0.3
        @test std(loodiff) * sqrt(n) ≈ 4.4 atol=0.2
    end
end


text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

const _rfft3_cache = Dict{Tuple{DataType,NTuple{3,Integer}},LinearOperator}()

function rfft3_operator(x::S)::LinearOperator where {S<:Array{T,3}} where {T<:AbstractFloat}
    key = (S, size(x))
    if !haskey(_rfft3_cache, key)
        y = copy(x)
        P = plan_rfft(y, flags = FFTW.MEASURE)
        forward(θ) = P * θ
        adjoint(θ) = P \ θ
        _rfft3_cache[key] = LinearOperator(size(x), size(x), forward, adjoint, false, T)
    end
    return _rfft3_cache[key]
end

function rfft3_operator(
    x::S,
)::LinearOperator where {S<:CuArray{T,3}} where {T<:AbstractFloat}
    key = (S, size(x))
    if !haskey(_rfft3_cache, key)
        y = copy(x)
        P = plan_rfft(y)
        forward(θ) = P * θ
        adjoint(θ) = P \ θ
        _rfft3_cache[key] = LinearOperator(size(x), size(x), forward, adjoint, false, T)
    end
    return _rfft3_cache[key]
end


text:
@testset "PredictiveAsymmetry" begin
    x, y, z = rand(100), rand(100), rand(100)
    ηs = 1:3

    @testset "Binning based" begin

        @testset "VisitationFrequency" begin
            # Define prediction lags and estimation method
            method = VisitationFrequency(RectangularBinning(3))

            # 𝔸(x → y) and  𝔸(x → y | z)
            𝔸reg  = predictive_asymmetry(x, y, method, ηs, normalize = false)
            𝔸cond = predictive_asymmetry(x, y, z, method, ηs, normalize = false)
            𝒜reg = predictive_asymmetry(x, y, method, ηs,  f = 1.0) # normalize == true by default
            𝒜cond = predictive_asymmetry(x, y, z, method, ηs, f = 1.5) # normalize == true by default

            @test 𝔸reg isa Vector{<:Real}
            @test 𝔸cond isa Vector{<:Real}
            @test 𝒜reg isa Vector{<:Real}
            @test 𝒜cond isa Vector{<:Real}
        end

        @testset "TransferOperator" begin
            # Define prediction lags and estimation method
            method = TransferOperator(RectangularBinning(2))

            # 𝔸(x → y) and  𝔸(x → y | z)
            𝔸reg  = predictive_asymmetry(x, y, method, ηs, normalize = false)
            𝔸cond = predictive_asymmetry(x, y, z, method, ηs, normalize = false)
            𝒜reg = predictive_asymmetry(x, y, method, ηs,  f = 1.0) # normalize == true by default
            𝒜cond = predictive_asymmetry(x, y, z, method, ηs, f = 1.5) # normalize == true by default

            @test 𝔸reg isa Vector{<:Real}
            @test 𝔸cond isa Vector{<:Real}
            @test 𝒜reg isa Vector{<:Real}
            @test 𝒜cond isa Vector{<:Real}
        end
    end

    @testset "Nearest neighbor based" begin
        @testset "Kraskov" begin
            # Define prediction lags and estimation method
            method = Kraskov()

            # 𝔸(x → y) and  𝔸(x → y | z)
            𝔸reg  = predictive_asymmetry(x, y, method, ηs, normalize = false)
            𝔸cond = predictive_asymmetry(x, y, z, method, ηs, normalize = false)
            𝒜reg = predictive_asymmetry(x, y, method, ηs,  f = 1.0) # normalize == true by default
            𝒜cond = predictive_asymmetry(x, y, z, method, ηs, f = 1.5) # normalize == true by default

            @test 𝔸reg isa Vector{<:Real}
            @test 𝔸cond isa Vector{<:Real}
            @test 𝒜reg isa Vector{<:Real}
            @test 𝒜cond isa Vector{<:Real}
        end

        @testset "KozachenkoLeonenko" begin
            # Define prediction lags and estimation method
            method = KozachenkoLeonenko()

            # 𝔸(x → y) and  𝔸(x → y | z)
            𝔸reg  = predictive_asymmetry(x, y, method, ηs, normalize = false)
            𝔸cond = predictive_asymmetry(x, y, z, method, ηs, normalize = false)
            𝒜reg = predictive_asymmetry(x, y, method, ηs,  f = 1.0) # normalize == true by default
            𝒜cond = predictive_asymmetry(x, y, z, method, ηs, f = 1.5) # normalize == true by default

            @test 𝔸reg isa Vector{<:Real}
            @test 𝔸cond isa Vector{<:Real}
            @test 𝒜reg isa Vector{<:Real}
            @test 𝒜cond isa Vector{<:Real}
        end
    end

    @testset "Kernel density based" begin
        @testset "NaiveKernel" begin
            # Define prediction lags and estimation method
            method = NaiveKernel(0.2)

            # 𝔸(x → y) and  𝔸(x → y | z)
            𝔸reg  = predictive_asymmetry(x, y, method, ηs, normalize = false)
            𝔸cond = predictive_asymmetry(x, y, z, method, ηs, normalize = false)
            𝒜reg = predictive_asymmetry(x, y, method, ηs,  f = 1.0) # normalize == true by default
            𝒜cond = predictive_asymmetry(x, y, z, method, ηs, f = 1.5) # normalize == true by default

            @test 𝔸reg isa Vector{<:Real}
            @test 𝔸cond isa Vector{<:Real}
            @test 𝒜reg isa Vector{<:Real}
            @test 𝒜cond isa Vector{<:Real}
        end
    end

    @testset "Hilbert" begin
        # Define prediction lags and estimation method
        method = Hilbert(VisitationFrequency(RectangularBinning(3)))

        # 𝔸(x → y) and  𝔸(x → y | z)
        𝔸reg  = predictive_asymmetry(x, y, method, ηs, normalize = false)
        𝔸cond = predictive_asymmetry(x, y, z, method, ηs, normalize = false)
        𝒜reg = predictive_asymmetry(x, y, method, ηs,  f = 1.0) # normalize == true by default
        𝒜cond = predictive_asymmetry(x, y, z, method, ηs, f = 1.5) # normalize == true by default

        @test 𝔸reg isa Vector{<:Real}
        @test 𝔸cond isa Vector{<:Real}
        @test 𝒜reg isa Vector{<:Real}
        @test 𝒜cond isa Vector{<:Real}
    end

    @testset "Symbolic" begin
        @testset "SymbolicPermutation" begin
            # Define prediction lags and estimation method
            method = SymbolicPermutation(m = 3)

            # 𝔸(x → y) and  𝔸(x → y | z)
            𝔸reg  = predictive_asymmetry(x, y, method, ηs, normalize = false)
            𝔸cond = predictive_asymmetry(x, y, z, method, ηs, normalize = false)
            𝒜reg = predictive_asymmetry(x, y, method, ηs,  f = 1.0) # normalize == true by default
            𝒜cond = predictive_asymmetry(x, y, z, method, ηs, f = 1.5) # normalize == true by default

            @test 𝔸reg isa Vector{<:Real}
            @test 𝔸cond isa Vector{<:Real}
            @test 𝒜reg isa Vector{<:Real}
            @test 𝒜cond isa Vector{<:Real}
        end
    end

end


text:
include(Pkg.dir("BismarkSummary","src","bismark-report.jl"))
test_bismark_report1 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report.txt")
test_bismark_report2 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report1.txt")
test_bismark_report3 = joinpath(Pkg.dir(), "BismarkSummary","testdata","bismark_report2.txt")
bismark_report_filenames = [test_bismark_report1,test_bismark_report2,test_bismark_report3]
sampleinfo_path = joinpath(Pkg.dir(), "BismarkSummary","testdata","datasource.tsv")
report_dict=parse_bismark_reports( bismark_report_filenames)
append_report_info_to_sampleinfo!(sampleinfo,report_dict)
plots = plot_bismark_summary(sampleinfo,["Sequence pairs analysed in total","Mapping efficiency"],["sex"])
hstack(plots...)


text:
# test durationstring output at borders
@test ProgressMeter.durationstring(0.9) == "0:00:00"
@test ProgressMeter.durationstring(1.0) == "0:00:01"
@test ProgressMeter.durationstring(59.9) == "0:00:59"
@test ProgressMeter.durationstring(60.0) == "0:01:00"
@test ProgressMeter.durationstring(60*60 - 0.1) == "0:59:59"
@test ProgressMeter.durationstring(60*60) == "1:00:00"
@test ProgressMeter.durationstring(60*60*24 - 0.1) == "23:59:59"
@test ProgressMeter.durationstring(60*60*24) == "1 days, 0:00:00"
@test ProgressMeter.durationstring(60*60*24*10 - 0.1) == "9 days, 23:59:59"
@test ProgressMeter.durationstring(60*60*24*10) == "10.00 days"


text:
# Copied from IsURL.jl
# Source: https://github.com/sindresorhus/is-absolute-url (MIT license)
const windowsregex = r"^[a-zA-Z]:[\\]"
const urlregex = r"^[a-zA-Z][a-zA-Z\d+\-.]*:"

"""
    isurl(str)
Checks if the given string is an absolute URL.
# Examples
```julia-repl
julia> isurl("https://julialang.org")
true
julia> isurl("mailto:someone@example.com")
true
julia> isurl("/foo/bar")
false
```
"""
function isurl(str::AbstractString)
    return !occursin(windowsregex, str) && occursin(urlregex, str)
end

function isgithub(url::AbstractString)
    if startswith(url, "http") && startswith(split(url, "//")[2], "github.com")
        return true
    else
        return false
    end
end


text:
# This file is a part of SimilaritySearch.jl

export JaccardDistance, DiceDistance, IntersectionDissimilarity
import Distances: evaluate

"""
    JaccardDistance()

The Jaccard distance is defined as

```math
J(u, v) = \\frac{|u \\cap v|}{|u \\cup v|}
```
"""
struct JaccardDistance <: PreMetric end

"""
    DiceDistance()

The Dice distance is defined as

```math
D(u, v) = \\frac{2 |u \\cap v|}{|u| + |v|}
```
"""
struct DiceDistance <: PreMetric end

"""
    IntersectionDissimilarity()

The intersection dissimilarity uses the size of the intersection as a mesuare of similarity as follows:

```math
I(u, v) = 1 - \\frac{|u \\cap v|}{\\max \\{|u|, |v|\\}}
```
"""
struct IntersectionDissimilarity <: PreMetric end

"""
    union_intersection(a::T, b::T)

Computes both the size of the unions an the size the intersections of `a` and `b`;
specified as ordered sequences.
"""
function union_intersection(a::T, b::T) where {T <: AbstractVector}
    len_a::Int = length(a)
    len_b::Int = length(b)
    ia::Int = 1
    ib::Int = 1
    intersection_size::Int = 0
    c::Int = 0
    @inbounds while ia <= len_a && ib <= len_b
        c = cmp(a[ia], b[ib])
        if c == 0
            ia += 1
            ib += 1
            intersection_size += 1
        elseif c < 0
            ia += 1
        else
            ib += 1
        end
    end

    len_a + len_b - intersection_size, intersection_size
end

"""
    evaluate(::JaccardDistance, a, b)

Computes the Jaccard's distance of `a` and `b` both sets specified as
sorted vectors.
"""
function evaluate(::JaccardDistance, a, b)
    u, i = union_intersection(a, b)
    1.0 - i / u
end

"""
    evaluate(::DiceDistance, a, b)

Computes the Dice's distance of `a` and `b` both sets specified as
sorted vectors.
"""
function evaluate(::DiceDistance, a, b)
    u, i = union_intersection(a, b)
    1.0 - 2 * i / (length(a) + length(b))
end

"""
    evaluate(::IntersectionDissimilarity, a, b)

Uses the intersection as a distance function (non-metric)
"""
function evaluate(::IntersectionDissimilarity, a, b)
    u, i = union_intersection(a, b)

    return 1.0 - i / max(length(a), length(b))
end


text: include("../../../../examples/histogram_msd.jl")


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   AbstractBandstructure{path}
#   --> path is the path in k-space
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractBandstructure{
    P<:AbstractReciprocalPath{RP} where {RP},
    H<:AbstractHamiltonian{L,UC,HB} where {L,UC,HB}
} end

# export the type
export AbstractBandstructure

################################################################################
#
#	INTERFACING / ACCESSING BANDSTRUCTURES
#	(functions have to be overwritten by concrete types)
#
################################################################################

# obtaining the hamiltonian
function hamiltonian(
            bs :: BS
        ) :: H where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, BS <: AbstractBandstructure{P,H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'hamiltonian' for bandstructure type " * string(typeof(bs)))
end

# export the function
export hamiltonian

# obtaining the path
function path(
            bs :: BS
        ) :: P where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, BS <: AbstractBandstructure{P,H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'path' for bandstructure type " * string(typeof(bs)))
end

# export the function
export path

# obtaining the energy values
# bands for each segment
# bands[i] gives all bands of segment i
# bands[i][j] gives all energy values for band j of segment i
# bands[i][j][k] gives the energy value at kpoint index k of band j in segment i
function energies(
            bs :: BS
        ) :: Vector{<:Vector{<:Vector{<:Real}}} where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, BS <: AbstractBandstructure{P,H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'energies' for bandstructure type " * string(typeof(bs)))
end

# export the function
export energies

# recalculate the band structure (energy values)
function recalculate!(
            bs :: BS
            ;
            resolution :: Union{Int64, Vector{Int64}} = 100
        ) where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, BS <: AbstractBandstructure{P,H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'recalculate!' for bandstructure type " * string(typeof(bs)))
end

# export the function
export recalculate!

function saveBandstructure(
        bs :: BS,
        fn :: AbstractString,
        group :: AbstractString = "bandstructure"
        ;
        append :: Bool = false
    ) where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, BS<:AbstractBandstructure{P,H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'recalculate!' for bandstructure type " * string(BS))
end

function loadBandstructure(
        ::Type{BS},
        fn :: AbstractString,
        group :: AbstractString = "bandstructure"
    ) where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}, BS<:AbstractBandstructure{P,H}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'recalculate!' for bandstructure type " * string(BS))
end

export saveBandstructure, loadBandstructure


text: ϵ = 0.0000001


text:
using CSV, Statistics, StatsBase
data = CSV.read("../data/temperatures.csv")[:,4]

println("Sample Mean: ", mean(data))
println("Harmonic <= Geometric <= Arithmetic ",
	(harmmean(data), geomean(data), mean(data)))
println("Sample Variance: ",var(data))
println("Sample Standard Deviation: ",std(data))
println("Minimum: ", minimum(data))
println("Maximum: ", maximum(data))
println("Median: ", median(data))
println("95th percentile: ", percentile(data, 95))
println("0.95 quantile: ", quantile(data, 0.95))
println("Interquartile range: ", iqr(data),"\n")

summarystats(data)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module Cond
using MLStyle

export @cond

function cond(cases, source, mod::Module)
    @match cases begin
        quote
        $(cases...)
        end =>
        let default = Expr(:call, throw, "None of the branches have satisfied conditions, at $(string(source)).")
            foldr(cases, init = default) do case, last
                last_lnode = source
                @match case begin
                    ::LineNumberNode => begin
                        last_lnode = case
                        Expr(:block, case, last)
                    end
                    :(_ => $b) => b
                    :($a => $b) => Expr(:if, a, b, last)
                    _ => throw("Invalid syntax for conditional branches at $last_lnode.")
                end
            end
        end
        _ => begin
            msg = "Malformed ast template, the second arg should be a block with a series of pairs(`a => b`), at $(string(source))."
            throw(SyntaxError(msg))
        end
    end
end

"""
```julia
@cond begin
    x > 1 => true
    _ => false
end
```
"""
macro cond(cases)
    cond(cases, __source__, __module__) |> esc
end

end


text:
function cartesian_iter(state, A, res, Asize)
    for i in CartesianIndices(Asize)
        idx = GPUArrays.gpu_sub2ind(Asize, i.I)
        res[idx] = A[idx]
    end
    return
end

function clmap!(state, f, out, b)
    i = linear_index(state) # get the kernel index it gets scheduled on
    out[i] = f(b[i])
    return
end

function ntuple_test(state, result, ::Val{N}) where N
    result[1] = ntuple(Val(N)) do i
        Float32(i) * 77f0
    end
    return
end

function ntuple_closure(state, result, ::Val{N}, testval) where N
    result[1] = ntuple(Val(N)) do i
        Float32(i) * testval
    end
    return
end

function test_base(AT)
    @testset "base functionality" begin
        @testset "copyto!" begin
            x = fill(0f0, (10, 10))
            y = rand(Float32, (20, 10))
            a = AT(x)
            b = AT(y)
            r1 = CartesianIndices((1:7, 3:8))
            r2 = CartesianIndices((4:10, 3:8))
            copyto!(x, r1, y, r2)
            copyto!(a, r1, b, r2)
            @test x == Array(a)
            r2 = CartesianIndices((4:11, 3:8))
            @test_throws DimensionMismatch copyto!(a, r1, b, r2)

            r2 = CartesianIndices((4:10, 3:8))
            x2 = fill(0f0, (10, 10))
            copyto!(x2, r1, b, r2)
            @test x2 == x

            fill!(a, 0f0)
            copyto!(a, r1, y, r2)
            @test Array(a) == x

            x = fill(0f0, (10,))
            y = rand(Float32, (20,))
            a = AT(x)
            b = AT(y)
            r1 = CartesianIndices((1:7,))
            r2 = CartesianIndices((4:10,))
            copyto!(x, r1, y, r2)
            copyto!(a, r1, b, r2)
            @test x == Array(a)
            r2 = CartesianIndices((4:11,))
            @test_throws ArgumentError copyto!(a, r1, b, r2)

            x = fill(0f0, (10,))
            y = rand(Float32, (20,))
            a = AT(x)
            b = AT(y)
            r1 = (1:7,)
            r2 = (4:10,)
            copyto!(x, r1[1], y, r2[1])
            copyto!(a, r1, b, r2)
            @test x == Array(a)

            x = fill(0., (10,))
            y = fill(1, (10,))
            a = AT(x)
            b = AT(y)
            copyto!(a, b)
            @test Float64.(y) == Array(a)
        end

        @testset "vcat + hcat" begin
            @test compare(vcat, AT, fill(0f0, (10, 10)), rand(Float32, 20, 10))
            @test compare(hcat, AT, fill(0f0, (10, 10)), rand(Float32, 10, 10))

            @test compare(hcat, AT, rand(Float32, 3, 3), rand(Float32, 3, 3))
            @test compare(vcat, AT, rand(Float32, 3, 3), rand(Float32, 3, 3))
            @test compare((a,b) -> cat(a, b; dims=4), AT, rand(Float32, 3, 4), rand(Float32, 3, 4))
        end

        @testset "reinterpret" begin
            a = rand(ComplexF32, 22)
            A = AT(a)
            af0 = reinterpret(Float32, a)
            Af0 = reinterpret(Float32, A)
            @test Array(Af0) == af0
            a = rand(ComplexF32, 4, 4)
            A = AT(a)
            @test_throws ArgumentError reinterpret(Int8, A)

            a = rand(ComplexF32, 10 * 10)
            A = AT(a)
            af0 = reshape(reinterpret(Float32, vec(a)), (20, 10))
            Af0 = reshape(reinterpret(Float32, vec(A)), (20, 10))
            @test Array(Af0) == af0
        end

        @testset "ntuple test" begin
            result = AT(Vector{NTuple{3, Float32}}(undef, 1))
            gpu_call(ntuple_test, result, (result, Val(3)))
            @test Array(result)[1] == (77, 2*77, 3*77)
            x = 88f0
            gpu_call(ntuple_closure, result, (result, Val(3), x))
            @test Array(result)[1] == (x, 2*x, 3*x)
        end

        @testset "cartesian iteration" begin
            Ac = rand(Float32, 32, 32)
            A = AT(Ac)
            result = fill!(copy(A), 0.0)
            gpu_call(cartesian_iter, result, (A, result, size(A)))
            Array(result) == Ac
        end

        @testset "Custom kernel from Julia function" begin
            x = AT(rand(Float32, 100))
            y = AT(rand(Float32, 100))
            gpu_call(clmap!, x, (-, x, y))
            jy = Array(y)
            @test map!(-, jy, jy) ≈ Array(x)
        end

        @testset "map" begin
            @test compare((a, b)-> map(+, a, b),    AT, rand(Float32, 10), rand(Float32, 10))
            @test compare((a, b)-> map!(-, a, b),   AT, rand(Float32, 10), rand(Float32, 10))
            @test compare((a, b, c, d)-> map!(*, a, b, c, d), AT, rand(Float32, 10), rand(Float32, 10), rand(Float32, 10), rand(Float32, 10))
        end

        @testset "repeat" begin
            @test compare(a-> repeat(a, 5, 6),  AT, rand(Float32, 10))
            @test compare(a-> repeat(a, 5),     AT, rand(Float32, 10))
            @test compare(a-> repeat(a, 5),     AT, rand(Float32, 5, 4))
            @test compare(a-> repeat(a, 4, 3),  AT, rand(Float32, 10, 15))
        end

        @testset "heuristics" begin
            blocks, threads = thread_blocks_heuristic(0)
            @test blocks == (1,)
            @test threads == (1,)
        end
    end
end


text:
# This is an image comparison via simple subtraction.
# All we do here is take the difference of each
# modeled trace with the real data trace. Whichever
# has the lowest error (and I'm using the term loosely,
# will be chosen by the FWI algorithm.
#
# author: GRAM 5/1/16

function compare(s_p, s_n, s_o, sei, i, j)

	# difference and weight
	dif_p = abs(sei[:,i,j] - s_p)
	dif_n = abs(sei[:,i,j] - s_n)
	dif_o = abs(sei[:,i,j] - s_o)

#=
	print("\n",dif_p,"\n")
	print("\n",dif_n,"\n")
	print("\n",dif_o,"\n")
=#

	# sum the difference vectors
	return [sum(dif_p),sum(dif_n),sum(dif_o)]

end


text:
using SimpleQuaternions, Test

a = SimpleQuaternion(1)
@test a==1

a = SimpleQuaternion(1,1,1,1)
@test a == 1 + im + jm + km
@test norm(a) == 2
@test a * inv(a) == 1
@test a/a == 1
@test a+a == 2a
@test a-a == 0
@test a*a == SimpleQuaternion(-2,2,2,2)
@test get_parts(a) == (1,1,1,1)
@test a//a == 1
@test sqrt(real(a*a')) == norm(a)

b = SimpleQuaternion(5,0,-2,4)
A = real_matrix(a)
B = real_matrix(b)
@test a == SimpleQuaternion(A)
@test a' == SimpleQuaternion(A')

@test A*B == real_matrix(a*b)

A = complex_matrix(a)
B = complex_matrix(b)
@test A*B == complex_matrix(a*b)
@test b == SimpleQuaternion(B)
@test b' == SimpleQuaternion(B')


text:
## -----------  By quality / parental

## Fraction in each quality, conditional on entry, by parental
# Counts returned are totals by parental, but returned by [qual, parental]
function frac_qual_by_parental(ds :: DataSettings)
    # This is enrollment by [quality, parental]
    dataV = read_matrix_by_xy(raw_entry_qual_parental(ds));
    @assert check_float_array(dataV, 0.005, 1.0);
    @check sum(dataV) ≈ 1.0
    nr = size(dataV, 1);

    # Make conditional on entry (columns sum to 1)
    # Now we have the fraction in each parental who attends each college (cond on entry)
    dataV = dataV ./ sum(dataV, dims = 1);
    @assert all(isapprox.(sum(dataV, dims = 1), 1.0))

    # Count in each cell
    cnts = read_matrix_by_xy(raw_entry_qual_parental(ds; momentType = :count));
    # Count by parental
    cnts = sum(cnts, dims = 1);
    cnts = repeat(cnts, outer = (nr, 1));

    ses = ses_from_choice_probs(dataV, cnts);
    cnts = round.(Int, cnts);
    return dataV, ses, cnts
end

## Graduation rates (conditional on entry) by (quality, parental)
function frac_grad_qual_parental(ds :: DataSettings)
    load_fct =
        mt -> read_matrix_by_xy(ds, :fracGrad_gpM, mt);
    m, ses, cnts = choice_prob_from_xy(load_fct);
    # rf = raw_frac_grad_qual_parental(ds);
    # dataM = read_matrix_by_xy(data_file(rf));
    # @assert check_float_array(dataM, 0.0, 1.0);
    @assert size(m) == (n_colleges(ds), n_parental(ds))
    # dataM = load_frac_grad_qual_parental(ds);
    # zero out 2 year colleges
    m[two_year_colleges(ds), :] .= 0.0;
    cnts[two_year_colleges(ds), :] .= 0;
    return m, ses, cnts
end

# function load_frac_grad_qual_parental(ds :: DataSettings)
#     rf = raw_frac_grad_qual_parental();
#     dataM = CollegeStratData.read_matrix_by_xy(CollegeStratData.data_file(rf));
#     @assert check_float_array(dataM, 0.0, 1.0);
#     @assert size(dataM) == (n_colleges(ds), n_parental(ds))
#     return dataM
# end

# 4y colleges only
function time_to_grad_4y_qual_parental(ds :: DataSettings)
    load_fct =
        mt -> read_matrix_by_xy(ds, :timeToGrad4y_qpM, mt);
    m, ses, cnts = load_mean_ses_counts(load_fct);
    @assert check_float_array(m, 3.0, 7.0);
    @assert size(m) == (n_4year(ds), n_parental(ds))
    # zero out 2 year colleges
    # m[two_year_colleges(ds), :] .= 0.0;
    # cnts[two_year_colleges(ds), :] .= 0;
    return m, ses, cnts
end

# function load_time_to_grad_qual_parental(ds :: DataSettings)
#     rf = raw_time_to_grad_qual_parental();
#     dataM = CollegeStratData.read_matrix_by_xy(CollegeStratData.data_file(rf));

#     @assert check_float_array(dataM, 3.0, 7.0);
#     @assert size(dataM) == (n_colleges(ds), n_parental(ds))
#     return dataM
# end

# ---------------


text:
#=
Example run from command line reading from the parameter file  examples/is_example?.jl (where ? = 1, 2, 3):
>  julia run.jl examples/is_example1

This run produces the output CSV file:  examples/in_example?.csv

This file contains functions for the infinite sites model in support of the Nearly Neutral paper
  ??? by Cheyenne Laue, ???, and Alden Wright
Author:  Alden H. Wright, Department of Computer Science, Univeristy of Montana, Missoula, MT 59812 USA
    alden.wright@umontana.edu
=#
using Dates
include("../src/InfSites.jl")

mutable struct infs_result_type
  nn_simtype::Int64
  N::Int64
  N_mu::Float64
  ngens::Int64
  dfe::Function
  dfe_str::AbstractString
  number_mutations::Int64
  number_extinctions::Int64
  number_fixations::Int64
  average_time_to_extinction::Float64
  average_time_to_fixation::Float64
  average_fitness_extinct::Float64
  average_fitness_fixed::Float64
  average_fitness_all::Float64
  stderr_sites_per_gen::Float64
  average_sites_per_gen::Float64
  average_heterozygosity_per_gen::Float64
  stderr_heterozygosity_per_gen::Float64
  count_fixed_del::Int64
  count_fixed_adv::Int64
end

function infs_result( nn_simtype::Int64, N::Int64, N_mu::Float64, ngens::Int64,
     dfe::Function=dfe_neutral, dfe_str::AbstractString="neutral" )
  tr = infs_result_type( nn_simtype, N, N_mu, ngens, dfe, dfe_str, 0, 0, 0, 0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0,0 )
  tr
end

function print_infs_result( tr::infs_result_type )
  if tr.nn_simtype == 1
    println("infinite alleles model")
  elseif tr.nn_simtype == 2
    println("infinite sites model")
  else
    println("nn_simtype: ", tr.nn_simtype)
  end
  println("N: ", tr.N)
  if mu_list_flag
    println("mu: ", tr.N_mu/tr.N)
  else
    println("N_mu: ", tr.N_mu)
  end
  println("ngens: ", tr.ngens)
  println("dfe: ", tr.dfe)
  println("dfe_str: ", tr.dfe_str)
  println("number_extinctions: ", tr.number_extinctions)
  println("number_fixations: ", tr.number_fixations)
  println("count_del_fixed: ", tr.count_fixed_del)
  println("count_adv_fixed: ", tr.count_fixed_adv)
  println("average time to extinction: ", tr.average_time_to_extinction)
  println("average time to fixation: ", tr.average_time_to_fixation)
  println("average fitness_extinct: ", tr.average_fitness_extinct)
  println("average fitness_fixed: ", tr.average_fitness_fixed)
  println("average fitness_all: ", tr.average_fitness_all)
  println("avg per generation count: ",tr.average_sites_per_gen)
  println("stderr per generation count: ",tr.stderr_sites_per_gen)
  println("avg per generation heterozygosity: ",tr.average_heterozygosity_per_gen)
  println("stderr per generation heterozygosity: ",tr.stderr_heterozygosity_per_gen)
  println()
end

function is_run_trials( mu_list_flag::Bool=false)
  println("stream: ",stream)
  trial = 1
  N = N_list[1]
  if !mu_list_flag
    tr = infs_result( nn_simtype, N, N_mu_list[1], ngens, dfe, dfe_str )
    is_writeheader(Base.stdout, N_list, tr, mu_list_flag=mu_list_flag )
    is_writeheader(stream, N_list, tr, mu_list_flag=mu_list_flag )
    for N_mu in N_mu_list
      for N in N_list
        tr = infs_result( nn_simtype, N, N_mu, ngens, dfe, dfe_str )
        is_run_trial( tr )
        is_writerow(Base.stdout, trial, tr )
        is_writerow(stream, trial, tr )
        Base.flush(stream)
        #print_infs_result( tr )
        Base.flush(Base.stdout)
        trial += 1
      end
    end
  else  # if mu_list_flag
    tr = infs_result( nn_simtype, N, N*mu_list[1], ngens, dfe, dfe_str )
    is_writeheader(Base.stdout, N_list, tr, mu_list_flag=mu_list_flag )
    is_writeheader(stream, N_list, tr, mu_list_flag=mu_list_flag )
    for mu in mu_list
      for N in N_list
        N_mu = N*mu
        tr = infs_result( nn_simtype, N, N_mu, ngens, dfe, dfe_str )
        is_run_trial( tr )
        is_writerow(Base.stdout, trial, tr, mu_list_flag=mu_list_flag )
        is_writerow(stream, trial, tr, mu_list_flag=mu_list_flag )
        Base.flush(stream)
        #print_infs_result( tr )
        Base.flush(Base.stdout)
        trial += 1
      end
    end
end  # if !mu_list_flag
end

function is_run_trial( tr::infs_result_type )
  if tr.nn_simtype == 2
    sc = InfSites.inf_sites( tr.N, tr.N_mu, tr.ngens, dfe=tr.dfe )
    tr.number_extinctions = length(sc.extinct)
    tr.number_fixations = length(sc.fixed)
    tr.count_fixed_del = sc.count_fixed_del
    tr.count_fixed_adv = sc.count_fixed_adv
    tr.average_time_to_extinction = InfSites.average_time_to_extinction(sc)
    tr.average_time_to_fixation = InfSites.average_time_to_fixation(sc)
    tr.average_fitness_extinct = InfSites.average_fitness_extinct(sc)
    tr.average_fitness_fixed = InfSites.average_fitness_fixed(sc)
    tr.average_fitness_all = InfSites.average_fitness_all(sc)
    (tr.average_sites_per_gen, tr.stderr_sites_per_gen) = InfSites.sites_per_gen(sc)
    (tr.average_heterozygosity_per_gen, tr.stderr_heterozygosity_per_gen) = InfSites.heterozygosity_per_gen(sc)
    return tr
  else
    println("nn_simtype ",nn_simtype," not implemented")
  end
end

@doc """ function writeheader(), slat_reps::Int64=100000 )
Write header line for output CSV file that corresponds to stream.
See the comments for run_simulation for a description of the parameters.
"""
function is_writeheader(stream::IO, N_list::Vector{Int64},
    tr::infs_result_type; mu_list_flag::Bool=false )
  dfe_params = [:dfe_adv_prob, :dfe_adv_alpha, :dfe_adv_theta, :dfe_disadv_prob, :dfe_disadv_alpha, :dfe_disadv_theta]
  param_strings = [
    "# $(string(Dates.today()))",
    (isdefined(Main,:seed) ? "# random number seed: $(seed)" : "# no random number seed defined"),
    #"# $((tr.nn_simtype==1) ? "infinite alleles model" :"infinite_sites_model")",
    "# infinite_sites_model",
    "# N_list=\"$(N_list)\"",
    (mu_list_flag ? "# mu_list=\"$(mu_list)\"" : "# N_mu_list=\"$(N_mu_list)\""),
    "# ngens=$(tr.ngens)",
    "# dfe=$(tr.dfe)",
    "# dfe_str=$(tr.dfe_str)"
    ]
  for dp in dfe_params
    if isdefined(Main,dp)
      Base.push!(param_strings,"# $(string(dp))=$(eval(dp))")
    end
  end
  write(stream, join(param_strings, "\n"), "\n")
  first_heads = [
    "type",
    #  "n",
    "N", (mu_list_flag ? "mu" : "N_mu")]
  mid_heads = []
  last_heads =
  [ "num_extinct",
    "num_fixed",
    "num_adv_fixed",
    "num_del_fixed",
    "fraction_fixed",
    "ave_extinct_time",
    "ave_fixed_time",
    "ave_extinct_selcoef",
    "ave_fixed_selcoef",
    "ave_all_selcoef",
    "ave_innov_count",
    "stderr_innov_count",
    "ave_heteroz",
    "stderr_heteroz",
  ]
  line = join(vcat( first_heads, mid_heads, last_heads), ",")
  write(stream, line, "\n")
end

function is_writerow(stream::IO, trial::Int64, tr::infs_result_type; mu_list_flag::Bool=false )
  first = Any[
    isdefined(Main,:type_str0) ? type_str : trial,
    #tr.n,           # sample size
    tr.N,           # popsize
    mu_list_flag ? tr.N_mu/tr.N : tr.N_mu,
  ]
 if tr.nn_simtype == 2
    mid = Any[
      tr.number_extinctions,
      tr.number_fixations,
      tr.count_fixed_adv,
      tr.count_fixed_del,
      Float64(tr.number_fixations)/(tr.number_extinctions+tr.number_fixations),
      tr.average_time_to_extinction,
      tr.average_time_to_fixation,
      tr.average_fitness_extinct,
      tr.average_fitness_fixed,
      tr.average_fitness_all,
      tr.average_sites_per_gen,
      tr.stderr_sites_per_gen,
      tr.average_heterozygosity_per_gen,
      tr.stderr_heterozygosity_per_gen
    ]
  else
    println("Error:  nn_simtype must be 2 for inf sites.")
    mid = Any[]
  end
  line = join( vcat( first, mid ), "," )
  write(stream, line, "\n")
end

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using VarianceComponentsHDFE
using Random, BenchmarkTools, Test
using JLD, SparseArrays, LinearAlgebra
# Make sure LDLFactorizations is version 0.5.0 and that the `multiple-rhs` branch is checked out

pkg_dir = pkgdir(VarianceComponentsHDFE)

# Obtain environment variables
run_large_benchmark = get(ENV, "VCHDFE_LARGE_BENCHMARK", "false")  == "true" ? true : false
num_rhs_large = try
        parse(Int, get(ENV, "BENCHMARK_NUM_RHS_LARGE", "2"))
    catch
        2
    end

# Entry point for the PkgBenchmarks call.  Can split into different files later.
include("prepare_benchmark_data.jl")
# NOTE: Suite below can assume that the `benchmark/data/...` has been filled

Random.seed!(1234)

# 1) Medium sized network benchmarks

medium_data = load(pkg_dir*"/benchmark/data/medium_main.jld")

X = medium_data["X_GroundedLaplacian"]
S_xx = medium_data["S_xx"]
S_xx_sparse = sparse(S_xx)

#Loading Operator
ldli, la = computeLDLinv(S_xx_sparse)
buffs = zeros(size(la)[1],1)
P = approxcholOperator(ldli, buffs[:,1])
compute_sol = approxcholSolver(P,la)

#RHS
m,k = size(X)
R_p = convert(Array{Float64,2}, bitrand(1,m))
rademacher!(R_p)
JLA_RHS = (R_p*X)[1,:]

#Benchmarks
@btime compute_sol($JLA_RHS;verbose=false)

#2) Large network benchmark
if run_large_benchmark

large_data = load(pkg_dir*"/benchmark/data/large_main.jld")

X = large_data["X_GroundedLaplacian"]
S_xx = large_data["S_xx"]
S_xx_sparse = sparse(S_xx)

#Loading Operator
ldli, la = computeLDLinv(S_xx_sparse)
buffs = zeros(size(la)[1],1)
P = approxcholOperator(ldli, buffs[:,1])
compute_sol = approxcholSolver(P,la)

#RHS
m,k = size(X)
R_p = convert(Array{Float64,2}, bitrand(1,m))
rademacher!(R_p)
JLA_RHS = (R_p*X)[1,:]

#Benchmarks
@btime compute_sol($JLA_RHS;verbose=false)

end


text:
## Simulation
using MonteCarlo, Distributions, DataFrames, JLD2

Tdist = Normal(MonteCarlo.IsingTc, .64)
n_Ts = 2^8
Ts = sort!(rand(Tdist, n_Ts))
Ts = Ts[Ts.>=1.2]
Ts = Ts[Ts.<=3.8]
therm = 10^4
sweeps = 10^3

df = DataFrame(L=Int[], T=Float64[], M=Float64[], χ=Float64[], E=Float64[], C_V=Float64[])

for L in 2 .^ [3, 4, 5, 6]
	println("L = ", L)
	for (i, T) in enumerate(Ts)
		println("\t T = ", T)
		beta = 1/T
		model = IsingModel(dims=2, L=L)
		mc = MC(model, beta=beta)
		run!(mc, sweeps=sweeps, thermalization=therm, verbose=false)
		meas = MonteCarlo.measurements(mc)[:ME]
		push!(df, [L, T, mean(meas[:Magn].m), mean(meas[:Magn].chi), mean(meas[:Energy].e), mean(meas[:Energy].C)])
	end
	flush(stdout)
end

sort!(df, [:L, :T])
@save "ising2d.jld" df

## Plots
using StatsPlots
@load "ising2d.jld" df

pyplot(size = (1024,690))
default(
	frame = :box,
	grid = :none,
	marker = :circle,
	markersize = 4.7,
)
p1 = @df df scatter(:T, :E, group = (L=:L,), xlab = "Temperature", ylab = "Energy", legend = :topleft)
p2 = @df df scatter(:T, :M, group = (L=:L,), xlab = "Temperature", ylab = "Magnetization",)
	x = range(1.2, stop=MonteCarlo.IsingTc, length=100)
	y = (1 .- sinh.(2.0 ./ (x)).^(-4)).^(1/8)
	plot!(p2, x, y, color = :black, linestyle = :dash, lw = 1.5, markeralpha = 0.0, markerstrokecolor = false, label = "exact")
	plot!(range(MonteCarlo.IsingTc, stop=3.8, length=100), zeros(100), color = :black, linestyle = :dash, lw = 1.5, markeralpha = 0.0, markerstrokecolor = false, label = nothing)
p3 = @df df scatter(:T, :C_V, group = (L=:L,), xlab = "Temperature", ylab = "Specific heat", )
	vline!(p3, [MonteCarlo.IsingTc], color = :black, linestyle = :dash, lw = 1.5, markeralpha = 0.0, markerstrokecolor = false, label = "\$ T_C \$")
p4 = @df df scatter(:T, :χ, group = (L=:L,), xlab = "Temperature", ylab = "Susceptibility",)
	vline!(p4, [MonteCarlo.IsingTc], color = :black, linestyle = :dash, lw = 1.5, markeralpha = 0.0, markerstrokecolor = false, label = "\$ T_C \$")
p = plot(p1,p2,p3,p4)
savefig(p, "ising2d.pdf")
savefig(p, "ising2d.svg")

## PyPlot
#
# using PyPlot
# grps = groupby(df, :L)
#
# # plot results together
# fig, ax = subplots(2,2, figsize=(12,8))
#
# for g in grps
# 	L = g.L[1]
# 	ax[1].plot(g.T, g.E, "o", markeredgecolor="black", label="L=$L")
# 	ax[2].plot(g.T, g.C_V, "o", markeredgecolor="black", label="L=$L")
# 	ax[3].plot(g.T, g.M, "o", markeredgecolor="black", label="L=$L")
# 	ax[4].plot(g.T, g.χ, "o", markeredgecolor="black", label="L=$L")
# end
# ax[1].legend(loc="best")
# ax[1].set_ylabel("Energy")
# ax[1].set_xlabel("Temperature")
#
# ax[2].set_ylabel("Specific heat")
# ax[2].set_xlabel("Temperature")
# ax[2].axvline(x=MonteCarlo.IsingTc, color="black", linestyle="dashed", label="\$ T_c \$")
# ax[2].legend(loc="best")
#
# ax[3].set_ylabel("Magnetization")
# ax[3].set_xlabel("Temperature")
# # ax[3].axvline(x=MonteCarlo.IsingTc, linewidth=2.0, color="black", label="\$ T_c \$")
# x = range(1.2, stop=MonteCarlo.IsingTc, length=100)
# y = (1 .- sinh.(2.0 ./ (x)).^(-4)).^(1/8)
# ax[3].plot(x,y, "k--", label="exact")
# ax[3].plot(range(MonteCarlo.IsingTc, stop=3.8, length=100), zeros(100), "k--")
# ax[3].legend(loc="best")
#
# ax[4].set_ylabel("Susceptibility χ")
# ax[4].set_xlabel("Temperature")
# ax[4].axvline(x=MonteCarlo.IsingTc, color="black", linestyle="dashed", label="\$ T_c \$")
# ax[4].legend(loc="best")
# tight_layout()
# savefig("ising2d.pdf")


text:
mutable struct WSumBatchedCache{Ta,Tb,Tc}
    σᵢₙ::Ta
    out::Tc
    δℒ::Tb
    valid::Bool
end

cache(l::WSum, arr_T, in_T, in_sz, batch_sz)  =
    WSumBatchedCache(similar(l.c, Complex{real(eltype(l.c))}, length(l.c),batch_sz),
              similar(l.c, Complex{real(eltype(l.c))}, 1, batch_sz),
              similar(l.c, Complex{real(eltype(l.c))}, 1, length(l.c), batch_sz),
              false)

batch_size(c::WSumBatchedCache) = size(c.out, 2)

function (l::WSum)(c::WSumBatchedCache, x)
    # dot product in temp cache
    c.σᵢₙ .= x .* l.c
    c.out .= sum_autobatch(c.σᵢₙ)

    # Store the input for backpropagation
    σ = copyto!(c.σᵢₙ, x)

    return c.out
end

function backprop(∇, l::WSum, c::WSumBatchedCache, δℒ::Number)
    # compute the derivative
    ∇.c .= c.σᵢₙ

    # Backpropagate
    c.δℒ .= δℒ .* transpose(l.c)
    return c.δℒ
end


text:
using WikiText, Test

@testset "WikiText" begin
    for wikitext in [WikiText2, WikiText103, WikiText2Raw, WikiText103Raw]
        for corpus in [wikitext, wikitext()]
            @test isfile(trainfile(corpus))
            @test isfile(validfile(corpus))
            @test isfile(testfile(corpus))
        end
    end
end


text:
####################################
# The AbstractBlockArray interface #
####################################

"""
    abstract AbstractBlockArray{T, N} <: AbstractArray{T, N}

The abstract type that represents a blocked array. Types that implement
the `AbstractBlockArray` interface should subtype from this type.

** Typealiases **

* `AbstractBlockMatrix{T}` -> `AbstractBlockArray{T, 2}`

* `AbstractBlockVector{T}` -> `AbstractBlockArray{T, 1}`

* `AbstractBlockVecOrMat{T}` -> `Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}`
"""
abstract AbstractBlockArray{T, N} <: AbstractArray{T, N}
typealias AbstractBlockMatrix{T} AbstractBlockArray{T, 2}
typealias AbstractBlockVector{T} AbstractBlockArray{T, 1}
typealias AbstractBlockVecOrMat{T} Union{AbstractBlockMatrix{T}, AbstractBlockVector{T}}

block2string(b, s) = string(join(map(string,b), '×'), "-blocked ", Base.dims2string(s))
Base.summary(a::AbstractBlockArray) = string(block2string(nblocks(a), size(a)), " ", typeof(a))
Base.similar{T}(block_array::AbstractBlockArray{T}) = similar(block_array, T)
Base.linearindexing{BA <: AbstractBlockArray}(::Type{BA}) = Base.LinearSlow()

"""
    nblocks(A, [dim...])

Returns a tuple containing the number of blocks in a block array.  Optionally you can specify
the dimension(s) you want the number of blocks for.

```jlcon
julia> A =  BlockArray(rand(5,4,6), [1,4], [1,2,1], [1,2,2,1]);

julia> nblocks(A)
(2,3,4)

julia> nblocks(A, 2)
3

julia> nblocks(A, 3, 2)
(4,3)
```
"""
nblocks(block_array::AbstractBlockArray, i::Int) = nblocks(block_array)[i]

function nblocks{N}(block_array::AbstractBlockArray, i::Vararg{Int, N})
    if N == 0
        throw(error("nblocks(A) not implemented"))
    end
    b = nblocks(block_array)
    return ntuple(k-> b[i[k]], Val{N})
end

"""
    blocksize(A, inds...)

Returns a tuple containing the size of the block at block index `inds...`.

```jlcon
julia> A = BlockArray(rand(5,4,6), [1,4], [1,2,1], [1,2,2,1]);

julia> blocksize(A,1,3,2)
(1,1,2)

julia> blocksize(A,2,1,3)
(4,1,2)
```
"""
function blocksize{T, N}(X, A::AbstractBlockArray{T,N}, ::Vararg{Int, N})
    throw(error("blocksize for ", typeof(A), "is not implemented"))
end

"""
    getblock(A, inds...)

Returns the block at blockindex `inds...`. An alternative syntax is `A[Block(inds...)].
Throws a `BlockBoundsError` if this block is out of bounds.

```jlcon
julia> A = BlockArray(rand(2,3), [1,1], [2,1])
2×2-blocked 2×3 BlockArrays.BlockArray{Float64,2,Array{Float64,2}}:
 0.190705  0.798036  │  0.471299
 --------------------┼-----------
 0.770005  0.845003  │  0.0315575

julia> getblock(A, 2, 1)
1×2 Array{Float64,2}:
 0.770005  0.845003

julia> A[Block(1,2)]
1×1 Array{Float64,2}:
 0.471299
```
"""
function getblock{T, N}(A::AbstractBlockArray{T,N}, ::Vararg{Int, N})
    throw("getblock for ", typeof(A), "is not implemented")
end

"""
    getindex{T, N}(A::AbstractBlockArray{T,N}, i...::Enum)

Returns the block at the blockindex defined by the values of the enums `i`

```jlcon
julia> A = PseudoBlockArray(zeros(2,3), [1,1], [2,1]);

julia> @enum vars u = 1 v = 2

julia> A[Block(1,2)] = [3.0];

ulia> A[u, v]
1×1 Array{Float64,2}:
 3.0
```
"""
@propagate_inbounds function Base.getindex{T, N}(A::AbstractBlockArray{T,N}, i::Vararg{Enum, N})
    getblock(A, map(Int, i)...)
end

@propagate_inbounds function Base.getindex{T, N}(A::AbstractBlockArray{T,N}, i::Enum)
    getblock(A, Int(i))
end

@propagate_inbounds function Base.getindex{T, N}(A::AbstractBlockArray{T,N}, i::Enum, j::Enum)
    getblock(A, Int(i), Int(j))
end

"""
    getblock!(X, A, inds...)

Stores the block at blockindex `inds` in `X` and returns it. Throws a `BlockBoundsError` if the
attempted assigned block is out of bounds.

```jlcon
julia> A = PseudoBlockArray(rand(2,3), [1,1], [2,1])
2×2-blocked 2×3 BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}:
 0.2062    0.0238446  │  0.0505515
 ---------------------┼-----------
 0.744768  0.225364   │  0.23028

julia> x = zeros(1,2);

julia> getblock!(x, A, 2, 1);

julia> x
1×2 Array{Float64,2}:
 0.744768  0.225364
```
"""
function getblock!{T, N}(X, A::AbstractBlockArray{T,N}, ::Vararg{Int, N})
    throw("getblock! for ", typeof(A), "is not implemented")
end

@propagate_inbounds function getblock!{T, N}(X, A::AbstractBlockArray{T,N}, i::Enum)
    getblock!(X, A, Int(i))
end

@propagate_inbounds function getblock!{T, N}(X, A::AbstractBlockArray{T,N}, i::Enum, j::Enum)
    getblock!(X, A, Int(i), Int(j))
end

@propagate_inbounds function getblock!{T, N}(X, A::AbstractBlockArray{T,N}, i::Vararg{Enum, N})
    getblock!(X, A, map(Int, i)...)
end

"""
    setblock!(A, v, inds...)

Stores the block `v` in the block at block index `inds` in `A`. An alternative syntax is `A[Block(inds...)] = v`.
Throws a `BlockBoundsError` if this block is out of bounds.

```jlcon
julia> A = PseudoBlockArray(zeros(2,3), [1,1], [2,1]);

julia> setblock!(A, [1 2], 1,1);

julia> A[Block(2,1)] = [3 4];

julia> A
2×2-blocked 2×3 BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}:
 1.0  2.0  │  0.0
 ----------┼-----
 3.0  4.0  │  0.0
```
"""
function setblock!{T, N}(A::AbstractBlockArray{T,N}, v, ::Vararg{Int, N})
    throw("setblock! for ", typeof(A), "is not implemented")
end

"""
    setindex!{T, N}(A::AbstractBlockArray{T,N}, v, i...::Enum)

Sets the block at the blockindex defined by the values of the enums `i`

```jlcon
julia> A = PseudoBlockArray(zeros(2,3), [1,1], [2,1]);

julia> @enum vars u = 1 v = 2

julia> A[u,v] = [3.0];

julia> A
2×2-blocked 2×3 BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}:
 0.0  0.0  │  3.0
 ----------┼-----
 0.0  0.0  │  0.0
```
"""
@propagate_inbounds function Base.setindex!{T, N}(A::AbstractBlockArray{T,N}, v, i::Vararg{Enum, N})
    setblock!(A, v, map(Int, i)...)
end

@propagate_inbounds function Base.setindex!{T, N}(A::AbstractBlockArray{T,N}, v, i::Enum)
    setblock!(A, v, Int(i))
end

@propagate_inbounds function Base.setindex!{T, N}(A::AbstractBlockArray{T,N}, v, i::Enum, j::Enum)
    setblock!(A, v, Int(i), Int(j))
end

"""
    BlockBoundsError([A],[inds...])

Thrown when a block indexing operation into a block array, `A`, tried to access an out-of-bounds block, `inds`.
"""
immutable BlockBoundsError <: Exception
    a::Any
    i::Any
    BlockBoundsError() = new()
    BlockBoundsError(a::AbstractBlockArray) = new(a)
    BlockBoundsError(a::AbstractBlockArray, i::ANY) = new(a,i)
end

function Base.showerror(io::IO, ex::BlockBoundsError)
    print(io, "BlockBoundsError")
    if isdefined(ex, :a)
        print(io, ": attempt to access ")
            print(io, summary(ex.a))
        if isdefined(ex, :i)
            print(io, " at block index [")
            print_joined(io, ex.i, ',')
            print(io, ']')
        end
    end
end

"""
    blockcheckbounds(A, inds...)

Throw a `BlockBoundsError` if the specified block indexes are not in bounds for the given block array. Subtypes of `AbstractBlockArray` should
specialize this method if they need to provide custom block bounds checking behaviors.

```jlcon
julia> A = BlockArray(rand(2,3), [1,1], [2,1]);

julia> blockcheckbounds(A, 3,2)
ERROR: BlockBoundsError: attempt to access 2×2-blocked 2×3 BlockArrays.BlockArray{Float64,2,Array{Float64,2}} at block index [3,2]
 in blockcheckbounds(::BlockArrays.BlockArray{Float64,2,Array{Float64,2}}, ::Int64, ::Int64)
 in eval(::Module, ::Any) at ./boot.jl:226
```
"""
function blockcheckbounds{T, N}(A::AbstractBlockArray{T, N}, i::Vararg{Int, N})
    n = nblocks(A)
    for (k, idx) in enumerate(i)
        if idx <= 0 || idx > n[k]
            throw(BlockBoundsError(A, i))
        end
    end
    return
end

"""
    full(A)

Returns the full array stored in `A`.
Full is here not used in the sense of sparse vs dense but in blocked vs unblocked.

```jlcon
julia> A = BlockArray(rand(2,3), [1,1], [2,1])
2×2-blocked 2×3 BlockArrays.BlockArray{Float64,2,Array{Float64,2}}:
 0.770528  0.396896  │  0.443308
 --------------------┼----------
 0.857069  0.403512  │  0.915934

julia> full(A)
2×3 Array{Float64,2}:
 0.770528  0.396896  0.443308
 0.857069  0.403512  0.915934
```
"""
function full(A::AbstractBlockArray) end

"""
    Block(inds...)

A `Block` is simply a wrapper around a set of indices or enums so that it can be used to dispatch on. By
indexing a `AbstractBlockArray` with a `Block` the a block at that block index will be returned instead of
a single element.

```jlcon
julia> A = BlockArray(rand(2,3), [1,1], [2,1])
2×2-blocked 2×3 BlockArrays.BlockArray{Float64,2,Array{Float64,2}}:
 0.190705  0.798036  │  0.471299
 --------------------┼-----------
 0.770005  0.845003  │  0.0315575

julia> A[Block(1,2)]
1×1 Array{Float64,2}:
 0.471299
```
"""
immutable Block{N, T}
    n::NTuple{N, T}
end

Block{N, T}(n::Vararg{T, N}) = Block{N, T}(n)

# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
@propagate_inbounds Base.getindex{T}(block_arr::AbstractBlockArray{T,1}, block::Block{1}) = getblock(block_arr, block.n[1])
@propagate_inbounds Base.getindex{T}(block_arr::AbstractBlockArray{T,2}, block::Block{2}) = getblock(block_arr, block.n[1], block.n[2])
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
@propagate_inbounds Base.getindex{T,N}(block_arr::AbstractBlockArray{T,N}, block::Block{N}) = getblock(block_arr, block.n...)

# vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #
@propagate_inbounds Base.setindex!{T}(block_arr::AbstractBlockArray{T,1}, v, block::Block{1}) =  setblock!(block_arr, v, block.n[1])
@propagate_inbounds Base.setindex!{T}(block_arr::AbstractBlockArray{T,2}, v, block::Block{2}) =  setblock!(block_arr, v, block.n[1], block.n[2])
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ #
@propagate_inbounds Base.setindex!{T,N}(block_arr::AbstractBlockArray{T,N}, v, block::Block{N}) =  setblock!(block_arr, v, block.n...)


text:
using ScratchNeuralNetwork
using Test

@testset "ScratchNeuralNetwork.jl" begin
    # Write your tests here.
end


text:
### A Pluto.jl notebook ###
# v0.16.4

using Markdown
using InteractiveUtils

# ╔═╡ 6d565956-314f-11ec-0005-a12d4cf1fb59
begin
	using Pkg; Pkg.activate("../")
	using nevanlinna_ac
	using DelimitedFiles
	using Plots
	gr(size=(600,500),
		xtickfontsize=13,
		ytickfontsize=13,
		xguidefontsize=16,
		yguidefontsize=16,
		legendfontsize=12,
		dpi=100,
		grid=(:y, :gray, :solid, 1, 0.4));
	using Printf
	using StatsFuns, SpecialFunctions
	"Packages"
end

# ╔═╡ 93711dc9-3b99-4ce3-8b39-570b4a833f58
η = 0.001

# ╔═╡ 25c4b6ca-c100-488e-8f73-e514f001aac8
dω = 4π/799

# ╔═╡ 1ed9dc6a-2971-4c61-aa3c-b3f91f7cb34b
md"""
### $g = 1.0, β = 20$
"""

# ╔═╡ eaa15cb3-8958-4967-b3b1-51f509b83b8a
md"""
### $g = 1.5, β = 20$
"""

# ╔═╡ 8f5723b4-5871-46b9-a706-495b1f34de77
md"""
### $g = 2.0, β = 20$
"""

# ╔═╡ 424f1e78-bbb1-4d55-8c14-a64d2f61ab0f
md"""
## check sum rule
"""

# ╔═╡ 2ed0a50b-2ab5-4bb7-9dbf-3af1f92c89c4
function chi_div_w(w::Vector, path::String, n::Int64)
	η = 0.001
	g = make_input(path)
	gs = MasubaraGF(n, g.GF[length(g.GF)-n+1:end])
	return [spectrum_density(i,η,gs) for i in w]
end

# ╔═╡ 64b77123-bab8-48ed-9291-409f6158c8ba
begin
	n =  5
	ω = [ i for i in range(-10,10,step=1.e-2)]

	grange = [1.0, 1.5, 2.0]
	brange = [10, 20, 30, 40]
	chi =Matrix{Vector}(undef,4,3)
	for i = 1:4, j = 1:3
		path = @sprintf "../data/gt/gt28_g_%.1f_b_%i.txt" grange[j] brange[i]
		chi[i,j] = chi_div_w(ω, path, n)
	end
end

# ╔═╡ cc64bdcc-a924-43b6-afb3-adc491eb1d3a
begin
	plot(ω, chi[1,1], label="β = 10, g = 1.0")
	plot!(ω, chi[2,1], label="β = 20, g = 1.0")
	plot!(ω, chi[3,1], label="β = 30, g = 1.0")
	plot!(ω, chi[4,1], label="β = 40, g = 1.0")
	plot!(xlim=(-1,1))
end

# ╔═╡ 3fcbf32e-f5ab-4fc0-95a7-f40696ea5de6
csum = 0.01 .* sum.(chi)

# ╔═╡ 1aad6d9d-7943-4742-aead-e946f3cb0e42
# check sum rule

# ╔═╡ 6569707a-8c7c-4692-8162-07796a8c65f7
md"""
## Compare $G(\tau)$
"""

# ╔═╡ ac03e64f-4857-4634-a86d-bc351aa84d75
begin
	dtau_1p0_1 = readdlm("../data/TFIsing-Li/gtau_1.0_b20.txt")
	dtau_1p0_2 = readdlm("../data/gt/gf_t8_g_1.0_b_20.txt")
	dtau_1p0_3 = readdlm("../data/gt/gf_t28_g_1.0_b_20.txt")
	gtau_1p0 = zeros(401,3)
	gtau_1p0[:,1] = dtau_1p0_1[:,1]
	gtau_1p0[:,2] = dtau_1p0_2[:,2] .- dtau_1p0_1[:,2]
	gtau_1p0[:,3] = dtau_1p0_3[:,2] .- dtau_1p0_1[:,2]
end

# ╔═╡ 97ca7dd6-ea05-4f18-a9b2-d2ccca9f423f
begin
	plot(gtau_1p0[:,1], gtau_1p0[:,2], label="χ=8")
	plot!(gtau_1p0[:,1], gtau_1p0[:,3], label="χ=2 * 8")
	plot!(xlabel="τ/β", ylabel="G(τ) loss", title="g = 1.0, β=20")
end

# ╔═╡ 04674261-d673-4ccc-a3f6-fe87e6332375
begin
	dtau_1p5_1 = readdlm("../data/TFIsing-Li/gtau_1.5_b20.txt")
	dtau_1p5_2 = readdlm("../data/gt/gf_t8_g_1.5_b_20.txt")
	dtau_1p5_3 = readdlm("../data/gt/gf_t28_g_1.5_b_20.txt")
	gtau_1p5 = zeros(401,3)
	gtau_1p5[:,1] = dtau_1p5_1[:,1]
	gtau_1p5[:,2] = dtau_1p5_2[:,2] .- dtau_1p5_1[:,2]
	gtau_1p5[:,3] = dtau_1p5_3[:,2] .- dtau_1p5_1[:,2]
end

# ╔═╡ 976e673c-a58e-4f09-be4a-f09d3512ef23
begin
	plot(gtau_1p5[:,1], gtau_1p5[:,2], label="χ=8")
	plot!(gtau_1p5[:,1], gtau_1p5[:,3], label="χ=2 * 8")
	plot!(xlabel="τ/β", ylabel="G(τ) loss", title="g = 1.5, β=20")
end

# ╔═╡ eb93aae6-1b7d-4055-b35e-d3b543248a1a
begin
	dtau_2p0_1 = readdlm("../data/TFIsing-Li/gtau_2.0_b20.txt")
	dtau_2p0_2 = readdlm("../data/gt/gf_t8_g_2.0_b_20.txt")
	dtau_2p0_3 = readdlm("../data/gt/gf_t28_g_2.0_b_20.txt")
	gtau_2p0 = zeros(401,3)
	gtau_2p0[:,1] = dtau_2p0_1[:,1]
	gtau_2p0[:,2] = dtau_2p0_2[:,2] .- dtau_2p0_1[:,2]
	gtau_2p0[:,3] = dtau_2p0_3[:,2] .- dtau_2p0_1[:,2]
end

# ╔═╡ 515b2659-7624-44a7-8cfb-108fde0bf5fc
begin
	plot(gtau_2p0[:,1], gtau_2p0[:,2], label="χ=8")
	plot!(gtau_2p0[:,1], gtau_2p0[:,3], label="χ=2 * 8")
	plot!(xlabel="τ/β", ylabel="G(τ) loss", title="g = 2.0, β=20")
end

# ╔═╡ 524af3b4-05eb-4269-bd7e-3de6f2113783
begin
	g = 1.5
	a,b = divrem(g, 1.0)
	a = Integer(a)
	b = Integer(b * 10)

	"""
	#import data
	p = @sprintf "../data/TFIsing-Li/origional/g_%ip%i.txt" a b
	op = @sprintf "../data/TFIsing-Li/gt_%.1f_b20.txt" g

	d = readdlm(p) ; len = length(d)/2 |> Integer
	t = [i for i in range(0,200,step=0.5)]

	#export data [t re.gt im.gt]
	open(op,"w") do file
		for i = 1:len
			writedlm(file,[t[i] d[2*i-1] d[2*i]])
		end
	end
	"""

	"""
	p = @sprintf "../data/TFIsing-Li/origional/g_%ip%i_tau_N_501.txt" a b
	op = @sprintf "../data/TFIsing-Li/gtau_%.1f_b20.txt" g

	#import data
	d = readdlm(p) ; len = length(d)
	tau = [i for i in range(0,20,step=0.05)]

	#export data [τ/β val]
	open(op,"w") do file
		for i = 1:len
			writedlm(file,[tau[i]/20 d[i]])
		end
	end
	"""

	"""
	nval = "A"
	p = @sprintf "../data/TFIsing-Li/origional/g_%ip%i_%s.txt" a b nval
	op = @sprintf "../data/TFIsing-Li/%sw_%.1f_b20.txt" nval g

	#import data
	d = readdlm(p) ; len = length(d)

	#make omega
	o = [-i for i in range(-2π,2π,step = 4π/len)]

	#export data [ω val]
	open(op,"w") do file
		for i = 1:len
			writedlm(file,[o[i] d[i]])
		end
	end
	"""

	"reexport data"
end

# ╔═╡ 42720925-c2f4-4103-a0a3-1cb51bc680f5
function f(w::Real)
	res = 1 - exp(-20*w)
	return res
end

# ╔═╡ 4147ddc7-bc31-47aa-a02e-eb193382b83c
function spectral_density(J::Real, ω::Real, β::Real; η::Float64=0.05,Γ::Real=1.)
    # 2Imχ(ω)
    if J/Γ == 0.
        return 2π*sinh(Γ*β)/cosh(Γ*β)*(delta(ω-2Γ,η) - delta(ω+2Γ,η))
    elseif J/Γ == 1.
        # Ref: PRX.4.031008(2014) A6
        g0 = 0.858714569; zc = J^(-1/4)
        T = 1/β
        up = zc * g0 * β^(3/4)
        down = 2^(1/4) * √π * gamma(1/8) * gamma(5/8)
        fac = up/down
        res = sinh(ω/(2*T)) * abs(gamma(1/8 - 1im*ω/(2π*T)))^2
        return 2*fac*res
    else
        @error "No exact results. J/Γ should be 1. or 0."
    end
end

# ╔═╡ 0b403313-abc1-4e91-a1f2-c745f167740f
begin
	n1 =  20
	#Zi-Long Li data
	A1_li= readdlm("../data/TFIsing-Li/Aw_1.0_b20.txt")
	omega = A1_li[:,1]
	A1w_li = A1_li[:,2] ./ omega

	#theory sachdev
	A1_theory = [spectral_density(1.0, i, 20) for i in omega]

	#cmpo data
	g1 = make_input("../data/gt/gt28_g_1.0_b_20.txt")
	g1s = MasubaraGF(10, g1.GF[length(g1.GF) - n1 + 1:end])
	A1_cmpo = [spectrum_density(w,η,g1s)*w for w in omega]
	A1w_cmpo = [spectrum_density(w,η,g1s) for w in omega]
end

# ╔═╡ 8060aa27-1de7-4b40-9665-15ce4760dee4
begin
	plot(omega, A1_theory, line=(2,:black), label="theory")
	plot!(omega, A1_li[:,2] ,line=(:dash,2), label="Li")
	plot!(omega, A1_cmpo, line=(:dash,2), label=@sprintf "χ=2*8,n=%i" n1)
	plot!(xlim=(0,5), ylim=(0,10))
	plot!(xlabel="ω", ylabel="A(ω)",title="g=1.0, β=20")
end

# ╔═╡ 1848b1c3-b52f-4b5a-9944-4d79e91a59a3
begin
	plot(omega, A1w_li,line=(2), label="Li")
	plot!(omega, A1w_cmpo, line=(:dash,2), label=@sprintf "χ=2*8,n=%i" n1)
	plot!(xlim=(-1,1))
	plot!(xlabel="ω", ylabel="A(ω)/ω",title="g=1.0, β=20")
end

# ╔═╡ 83ad4c7b-8174-4d63-bb3c-33aa4faa3edb
begin
	n1p5 =  5
	#Zi-Long Li data
	A1p5_li= readdlm("../data/TFIsing-Li/Aw_1.5_b20.txt")
	A1p5w_li = A1p5_li[:,2] ./ omega

	#cmpo data
	g1p5 = make_input("../data/gt/gt28_g_1.5_b_20.txt")
	g1p5s = MasubaraGF(n1p5, g1p5.GF[length(g1p5.GF) - n1p5 + 1:end])
	A1p5_cmpo = [spectrum_density(w,η,g1p5s)*w for w in omega]
	A1p5w_cmpo = [spectrum_density(w,η,g1p5s) for w in omega]
end

# ╔═╡ fcb57525-938a-4951-8e21-11ca03a3b318
begin
	plot(omega, A1p5w_li ,line=(2), label="Li")
	plot!(omega, A1p5w_cmpo, line=(:dash,2), label=@sprintf "χ=2*8,n=%i" n1p5)
	plot!([-1,1], seriestype="vline",line=(:dash, :black))
	plot!(xlabel="ω", ylabel="A(ω)/ω",title="g=1.5, β=20")
end

# ╔═╡ 987bd147-f917-47a2-ad04-2571c9e69dd7
begin
	n2 =  5
	#Zi-Long Li data
	A2_li= readdlm("../data/TFIsing-Li/Aw_2.0_b20.txt")
	A2w_li = A2_li[:,2] ./ omega

	#cmpo data
	g2 = make_input("../data/gt/gt28_g_2.0_b_20.txt")
	g2s = MasubaraGF(n2, g2.GF[length(g2.GF) - n2 + 1:end])
	A2_cmpo = [spectrum_density(w,η,g2s)*w for w in omega]
	A2w_cmpo = [spectrum_density(w,η,g2s) for w in omega]
end

# ╔═╡ cef09d6b-9fd8-47c8-9b6d-4def6c6d0a2a
begin
	plot(omega, A2w_li ,line=(2), label="Li")
	plot!(omega, A2w_cmpo, line=(:dash,2), label=@sprintf "χ=2*8,n=%i" n2)
	plot!([-2,2], seriestype="vline",line=(:dash, :black))
	plot!(xlabel="ω", ylabel="A(ω)/ω",title="g=1.5, β=20")
end

# ╔═╡ 73760863-b866-4f56-8115-8b03ab7c8192
pwd()

# ╔═╡ Cell order:
# ╟─93711dc9-3b99-4ce3-8b39-570b4a833f58
# ╟─25c4b6ca-c100-488e-8f73-e514f001aac8
# ╟─1ed9dc6a-2971-4c61-aa3c-b3f91f7cb34b
# ╟─0b403313-abc1-4e91-a1f2-c745f167740f
# ╟─8060aa27-1de7-4b40-9665-15ce4760dee4
# ╟─1848b1c3-b52f-4b5a-9944-4d79e91a59a3
# ╟─eaa15cb3-8958-4967-b3b1-51f509b83b8a
# ╟─83ad4c7b-8174-4d63-bb3c-33aa4faa3edb
# ╠═fcb57525-938a-4951-8e21-11ca03a3b318
# ╟─8f5723b4-5871-46b9-a706-495b1f34de77
# ╟─987bd147-f917-47a2-ad04-2571c9e69dd7
# ╟─cef09d6b-9fd8-47c8-9b6d-4def6c6d0a2a
# ╟─424f1e78-bbb1-4d55-8c14-a64d2f61ab0f
# ╟─2ed0a50b-2ab5-4bb7-9dbf-3af1f92c89c4
# ╟─64b77123-bab8-48ed-9291-409f6158c8ba
# ╟─cc64bdcc-a924-43b6-afb3-adc491eb1d3a
# ╠═3fcbf32e-f5ab-4fc0-95a7-f40696ea5de6
# ╠═1aad6d9d-7943-4742-aead-e946f3cb0e42
# ╟─6569707a-8c7c-4692-8162-07796a8c65f7
# ╟─ac03e64f-4857-4634-a86d-bc351aa84d75
# ╟─97ca7dd6-ea05-4f18-a9b2-d2ccca9f423f
# ╟─04674261-d673-4ccc-a3f6-fe87e6332375
# ╟─976e673c-a58e-4f09-be4a-f09d3512ef23
# ╟─eb93aae6-1b7d-4055-b35e-d3b543248a1a
# ╟─515b2659-7624-44a7-8cfb-108fde0bf5fc
# ╟─524af3b4-05eb-4269-bd7e-3de6f2113783
# ╟─42720925-c2f4-4103-a0a3-1cb51bc680f5
# ╟─4147ddc7-bc31-47aa-a02e-eb193382b83c
# ╟─73760863-b866-4f56-8115-8b03ab7c8192
# ╠═6d565956-314f-11ec-0005-a12d4cf1fb59


text:
using Documenter, SolverTraces

makedocs(;
    modules=[SolverTraces],
    format=Documenter.HTML(assets=["assets/style.css"]),
    pages=[
        "Home" => "index.md",
        "Solver Trace Columns" => "columns.md",
        "Colors" => "colors.md",
    ],
    repo="https://github.com/jagot/SolverTraces.jl/blob/{commit}{path}#L{line}",
    sitename="SolverTraces.jl",
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com>",
    doctest=false
)

deploydocs(;
    repo="github.com/jagot/SolverTraces.jl",
)


text:
import Base: one
@deprecate one(::Type{NumDenom{T}}) where T     oneunit(NumDenom{T})
@deprecate one(p::NumDenom)                     oneunit(p)

@deprecate ratio(mm::AbstractArray, args...)  ratio.(mm, args...)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenEXR_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenEXR")
JLLWrappers.@generate_main_file("OpenEXR", UUID("18a262bb-aa17-5467-a713-aee519bc75cb"))
end  # module OpenEXR_jll


text:
immutable FirmThreshold{T <: Union{Real, AbstractArray}} <: ProximableFunction
  lambda::T
  tau::T
  function FirmThreshold(lambda::T, tau::T )

	  if !(eltype(lambda) <: Real)
	          error("λ must be real")
	  end
	  if !(eltype(tau) <: Real)
	          error("τ must be real")
	  end
	  if any(lambda .< 0)
	          error("λ must be nonnegative")
	  end
	  if any(tau .< 0) || any(tau .> 1)
		  error("τ must be 0 < τ < 1")
	  end
	  if length(tau) != length(lambda)
	          error("λ and τ must have the same length")
	  end
	  new(lambda,tau)
  end
end

"""
  FirmThreshold(λ::Real=1.0,τ::Real=0.7)

  Returns the function `g(x) = sum( τ(λ|x_i|-x_i^2/2) for |x_i|<λ; τλ^2/2 elsewhere)`, for a real parameter `λ ⩾ 0` and `0 < τ < 1`.
"""
FirmThreshold{R <: Real}(lambda::R=1.0,tau::R=0.7) = FirmThreshold{R}(lambda,tau)

FirmThreshold{A <: AbstractArray}(lambda::A,tau::A) = FirmThreshold{A}(lambda,tau)

function (f::FirmThreshold{R}){R <: Real}(x::AbstractArray)
  fy = 0.
  for i in eachindex(x)
	  fy += abs(x[i]) < f.lambda ? f.tau*(f.lambda*abs(x[i]) -0.5*real(dot(x[i],x[i]))) : 0.5*f.tau*f.lambda^2
  end
  return fy
end

function (f::FirmThreshold{A}){A <: AbstractArray}(x::AbstractArray)
  fy = 0.
  for i in eachindex(x)
	  fy += abs(x[i]) < f.lambda[i] ? f.tau[i]*(f.lambda[i]*abs(x[i]) - 0.5*real(dot(x[i],x[i]))) : 0.5*f.tau[i]*f.lambda[i]^2
  end
  return fy
end

function prox!{A <: AbstractArray, R <: RealOrComplex}(f::FirmThreshold{A}, x::AbstractArray{R}, y::AbstractArray{R}, gamma::Real=1.0)
  fy = zero(R)
  for i in eachindex(x)
	  y[i] = abs(x[i]) < f.lambda[i]*f.tau[i] ? 0 : ( abs(x[i]) >= f.lambda[i] ? x[i] : (x[i]-f.lambda[i]*f.tau[i]*sign(x[i]))/(1-f.tau[i])  )
  end
  return f(y)
end

function prox!{T <: Real, R <: RealOrComplex}(f::FirmThreshold{T}, x::AbstractArray{R}, y::AbstractArray{R}, gamma::Real=1.0)
  fy = zero(R)
  for i in eachindex(x)
	  y[i] = abs(x[i]) < f.lambda*f.tau ? 0 : ( abs(x[i]) >= f.lambda ? x[i] : (x[i]-f.lambda*f.tau*sign(x[i]))/(1-f.tau)  )
  end
  return f(y)
end

fun_name(f::FirmThreshold) = "weighted Firm Threshold"
fun_dom(f::FirmThreshold) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::FirmThreshold) = "x ↦ sum(τ(λ|x_i|-x_i^2/2) for |x_i|<λ; τλ^2/2 elsewhere)"
fun_params{R <: Real}(f::FirmThreshold{R}) = "λ = $(f.lambda)"
fun_params{A <: AbstractArray}(f::FirmThreshold{A}) = string("λ = ", typeof(f.lambda), " of size ", size(f.lambda))

function prox_naive{A<: Real,T <: RealOrComplex}(f::FirmThreshold{A}, x::AbstractArray{T}, gamma::Real=1.0)
  y = similar(x)
  for i in eachindex(x)
	  if abs(x[i]) < f.lambda*f.tau
		  y[i] = 0
	  elseif abs(x[i]) >= f.lambda
		  y[i] = x[i]
	  else
		  y[i] = (x[i]-f.lambda*f.tau*sign(x[i]))/(1-f.tau)
	  end
  end
  fy = 0.
  for i in eachindex(x)
          if abs(y[i]) < f.lambda
        	  fy += f.tau*(f.lambda*abs(y[i])-0.5*real(dot(y[i],y[i])))
          else
        	  fy += 0.5*f.tau*f.lambda^2
          end
  end
  return y, fy
end

function prox_naive{A<: AbstractArray,T <: RealOrComplex}(f::FirmThreshold{A}, x::AbstractArray{T}, gamma::Real=1.0)
  y = similar(x)
  for i in eachindex(x)
	  if abs(x[i]) < f.lambda[i]*f.tau[i]
		  y[i] = 0
	  elseif abs(x[i]) >= f.lambda[i]
		  y[i] = x[i]
	  else
		  y[i] = (x[i]-f.lambda[i]*f.tau[i]*sign(x[i]))/(1-f.tau[i])
	  end
  end
  fy = 0.
  for i in eachindex(x)
	  if abs(y[i]) < f.lambda[i]
		  fy += f.tau[i]*(f.lambda[i]*abs(y[i])-0.5*real(dot(y[i],y[i])))
          else
		  fy += 0.5*f.tau[i]*f.lambda[i]^2
          end
  end
  return y, fy
end


text:
function telegrapher(eq, dofs, dofsneigh, flux, fluxneigh, dx, normalidx, normalsign, numericalflux)
    if normalidx != 1
        numericalflux .= 0.0
        return
    end

    maxeigenval_left = max_eigenval(eq, dofs, normalidx)
    maxeigenval_right = max_eigenval(eq, dofsneigh, normalidx)
    maxeigenval = max(maxeigenval_left, maxeigenval_right)

    s = TelegraphersShortcuts()
    C = dofs[s.C]
    L = dofs[s.L]
    Il = dofs[s.I]
    Ir = dofsneigh[s.I]
    Vl = dofs[s.V]
    Vr = dofsneigh[s.V]
    Istar = 0.5 .* (Il .+ Ir) + 0.5 * sqrt(C/L) * normalsign * (Vl - Vr)
    Vstar = 0.5 .* (Vl .+ Vr) + 0.5 * sqrt(L/C) * normalsign * (Il - Ir)

    numericalflux[s.I] = 1/L * Vstar
    numericalflux[s.V] = 1/C * Istar
end


text:
using Pkg
Pkg.API.develop(Pkg.PackageSpec(name="MrPhelps", path="/home/caseykneale/Desktop/Playground/MrPhelps/"))
using MrPhelps
using CSV, DataFrames, Random

const entries = 111

datadir  = "Playground/MrPhelps/data/"
filefmt  = "expr_{experiment}_number_{num}.csv"
fullpath = Base.joinpath(datadir, filefmt)

fauxfilenames = Expand( fullpath, Dict( "experiment" => [ "A", "B", "C" ],
                                        "num"     => [ "1", "2" ] )
                      )

#Make some random 3 column CSV's and save them to disk via the expand call
for filename in fauxfilenames
    column_names = [ :name, :rand, :randn ]
    column_values = [ [ randstring(rand(1:10)[1] + 4)   for entry in 1:entries ],
                      [ rand(1)[1]                      for entry in 1:entries ],
                      [ randn(1)[1]                     for entry in 1:entries ] ]
    df = DataFrame( column_values, column_names )
    CSV.write(filename, df)
end

mock_data_meta = VariableGlob( datadir, filefmt )

println.( mock_data_meta.parsed_values )

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps