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

using Markdown
using InteractiveUtils

# ╔═╡ feb1b280-268f-11ec-1574-e5037552b641
begin

# Import packages used

using Distributions
using DataFrames
using Optim
using CSV
using DelimitedFiles
using StatFiles
using Statistics
using HypothesisTests
using Distributions

end

# ╔═╡ d08aaa10-268f-11ec-067f-4322e82ba779
# Augenblick and Rabin, 2019, "An Experiment on Time Preference and Misprediction in Unpleasant Tasks", Table 2

# ╔═╡ fee9d980-268f-11ec-3901-3705290b7c24
# Authors:

# Massimiliano Pozzi (Bocconi University, pozzi.massimiliano@studbocconi.it)
# Salvatore Nunnari (Bocconi University, https://snunnari.github.io, salvatore.nunnari@unibocconi.it)

# The code in this Jupyter notebook performs the individual estimates to replicate columns 1, 2, 3, 4 of Table 2

# This notebook was tested with the following packages versions:

# Pozzi: julia 1.6.1, Pluto 0.14.7, Neptune 0.14.0, Distributions 0.25.17, DataFrames 1.2.2, Optim 1.4.1, CSV 0.9.5

# ╔═╡ f8c1b410-51b4-11ec-1d65-11df2949bf23
# This notebook works with both Neptune and Pluto; the "begin" and "end" tags in each cell are necessary for Pluto notebooks but do not play any role in Neptune

# ╔═╡ 0de27c80-51b5-11ec-05e9-a5b37469b0c7
# To import a new package (es. StatFiles) use this cell. Change the name inside the brackets to change Pkg

# import Pkg; Pkg.add("StatFiles")

# ╔═╡ fe9ddc60-268f-11ec-112a-59ea860b2c81
begin

# 1. Data Cleaning and Data Preparation

# We import the dataset containing observations for all 100 individuals. We then remove observations when a bonus was offered and create dummy variables that will be useful for estimation: pb is equal to one if the subject completed 10 mandatory tasks on subject-day (this is used to estimate the projection bias parameter α), ind_effort10 and ind_effort110 are equal to one if the subject completed 10 or 110 tasks respectively (and they are used for the Tobit correction when computing the likelihood)

dt = DataFrame(load(raw"../input/decisions_data.dta")) # full sample
dt = dt[dt.bonusoffered .!=1,:]  # remove observations when a bonus was offered

dt.pb = dt.workdone1 ./ 10 # pb dummy variable. workdone1 can either be 10 or 0, so                                dividing the variable by 10 creates our dummy

dt.ind_effort10  = dt.effort .== 10    # ind_effort10 dummy
dt.ind_effort110 = dt.effort .== 110   # ind_effort110 dummy

end;

# ╔═╡ 74270f10-51b5-11ec-2edd-af650c9bec6b
# 2. Define the Model and the Likelihood (Section 3 in the Paper)

# For a more detailed explanation of the model please refer to the python notebook

# Define the criterion function to minimize. This function computes the negative log likelihood of observing the choices made by the agents. We first compute the predicted choice that is found by solving the agents' utility maximization problem. We then assume that the observed effort is distributed as a Gaussian with mean the predicted choice and a given value of sigma (to be estimated). We then work out the likelihood of observing the effort found in the data, take the log and the negative of it to define our objective function.

# parameters:

# beta is the present bias parameter
# betahat is the perceived present bias parameter
# delta is the usual time-discounting parameter
# gamma and phi are the two parameters controlling the cost of effort function
# alpha is the projection bias parameter
# sigma is the standard deviation of the normal error term ϵ

# args:

# netdistance is the difference between the payment date T and the work time t
# wage is the amount paid per task in a certain session
# today is a dummy variable equal to one if the decision involves the choice of work today
# prediction is a dummy variable equal to one if the decision involves the choice of work in the future
# pb is a dummy equal to one if the subject completed 10 mandatory tasks on subject-day
# effort is the number of tasks completed by a subject in a session. It can range from a minimum of 10 to a maximum of 110
# ind_effort10 is a dummy equal to one if the subject's effort was equal to 10
# ind_effort110 is a dummy equal to one if the subject's effort was equal to 110

function negloglike(params, args)

    netdistance, wage, today, prediction, pb, effort, ind_effort10, ind_effort110, spec = args

# Column 4 specification has 7 parameters, the others have 6 parameters. Use if statement to specify the correct params vector

	if spec!=4
		beta, betahat, delta, gamma, phi, sigma  = params
	else
		beta, betahat, delta, gamma, phi, alpha, sigma = params
	end

	# compute the predicted choice. To avoid domain error we use max for beta,               betahat, delta, phi

	predchoice=((max(phi,1e-6).*(max(delta,1e-6).^netdistance).*(max(beta,1e-6).^today).*(max(betahat,1e-6).^prediction).*wage).^(1 ./(gamma.-1)))

	if spec==4 predchoice = predchoice .- pb.*alpha end

	# compute the vector of individual observations' likelihood, a 1x8049 vector

    prob = pdf.(Normal(0,max(sigma,0.1)),predchoice.-effort)

	# find index when effort is 10 ot 110 to apply Tobit correction
	len = length(effort)
    index_10  = [i for i=1:len if effort[i]==10]
    index_110 = [i for i=1:len if effort[i]==110]

	# Apply Tobit correction
    for i in index_10
        prob[i] = 1 - cdf.(Normal(0,1),(predchoice[i]-effort[i])/max(sigma,0.1))
	end

    for i in index_110
        prob[i] = cdf.(Normal(0,1),(predchoice[i]-effort[i])/max(sigma,0.1))
	end

	# find index when prob is equal to 0 or 1 and change it. this is needed to avoid         problems when taking logs
    index_p0 = [i for i=1:len if prob[i]==0]
    index_p1 = [i for i=1:len if prob[i]==1]

	# change prob when 0 or 1
    for i in index_p0
        prob[i] = 1e-6
	end

    for i in index_p1
        prob[i] = 1 - 1e-6
	end

    negll = - sum(log.(prob)) # negative log likelihood

    return negll

end

# ╔═╡ ce679cb0-51b5-11ec-39ba-3f65d7b1ce33
# Estimation

# We define the function that computes the individual estimates for columns 1, 2, 3 and 4 of Table 2. There are four specifications to consider: the first column uses the entire sample, the second column uses only observations for which the decision date is less than 4, the third column uses only observations for which the decision date is from time 4 onwards, and the fourth column uses the full sample but estimates also the projection bias parameter.

# Define an auxiliary function.
# This function takes as input a dataset, a scalar and a specification argument and returns the relevant individual-level dataframe containing only observations for individual whose ID is equal to that scalar for a given specification

function getind(dataset, wid, spec)
    dataset_ind = dataset[dataset.wid .== wid,:]
	# second column specification
	if spec == 2 dataset_ind = dataset_ind[dataset_ind.decisiondatenum .< 4,:]  end
	# third column specification
	if spec == 3 dataset_ind = dataset_ind[dataset_ind.decisiondatenum .>= 4,:] end
    return dataset_ind

end

# ╔═╡ a5d0925e-5535-11ec-06ca-85bfd5b5081b
# Define the function that estimates the model

function estimate(dataset, spec)

# Define the vectors containing the initial guesses and the vectors that will contain our results. res_params will be a 1x100 vector, each element is a 1x6 or 1x7 vector containing the estimates for the parameters. res_like will be a 1x100 vector containing the log likelihood obtained, while res_wid will be a 1x100 vector containing the wid of the individual for which we find the estimates. res_conv is a 1x100 vector of 1 (converged) or 0 (not converges)

    beta_init, betahat_init, delta_init, gamma_init, phi_init, alpha_init, sigma_init = 1.0, 1, 1, 2, 250, 3, 50

	res_param = []
    res_like  = Vector{Float64}(undef, 100)
    res_wid   = Vector{Float64}(undef, 100)
	res_conv  = Vector{Float64}(undef, 100)

# compute the individual estimates and save the results

    for wid in unique(dataset.wid)  # loop over individuals' ID

        params_init = [beta_init, betahat_init, delta_init, gamma_init, phi_init,                            alpha_init, sigma_init]
        params_nopb = [beta_init, betahat_init, delta_init, gamma_init, phi_init,                            sigma_init]

        dt_ind = getind(dataset, wid, spec) # get the relevant individual                                                           dataframe that depends also on the                                                     specification

# args needed to compute the individual likelihood

        args_ind=[dt_ind.netdistance, dt_ind.wage, dt_ind.today, dt_ind.prediction, dt_ind.pb, dt_ind.effort, dt_ind.ind_effort10, dt_ind.ind_effort110, spec]

# We run each estimation twice starting from the previous argmin to see if we were stuck on a local minima. This is usually not necessary, but also not expensive in terms of time so it can be worth it to do it

        if spec != 4 # first three columns do not have the projection bias parameter

            j = 0
            while (j<3)
                global sol = optimize(x->negloglike(x,args_ind),
				   params_nopb,
				   method = NelderMead(),
				   iterations = 2500)
				   params_nopb = Optim.minimizer(sol)
                j = j+1
			end
        else
            j = 0
            while (j<3)
                global sol = optimize(x->negloglike(x,args_ind),
				   params_init,
				   method = NelderMead(),
				   iterations = 4000) # projetion bias needs more iterations
				   params_init = Optim.minimizer(sol)
                j = j+1
			end
		end

    append!(res_param, [Optim.minimizer(sol)])
    res_like[round(Int32, wid)]  = -Optim.minimum(sol)
    res_wid[round(Int32, wid)]   = wid
	res_conv[round(Int32, wid)]  = Optim.converged(sol)

	end # end for wid loop

    return res_param, res_like, res_conv

end

# ╔═╡ f2177020-5536-11ec-0db6-1b2a59aa1c0e
begin

# This cell computes the estimates for column (1) and stores the results.

res1 = estimate(dt,1)
param1 = res1[1]
like1  = res1[2]
conv1  = res1[3]

end;

# ╔═╡ 99348ad0-5543-11ec-274d-1d59ea3e02b0
begin

# This cell computes the estimates for column (2) and stores the results.

res2 = estimate(dt,2)
param2 = res2[1]
like2  = res2[2]
conv2  = res2[3]

end;

# ╔═╡ a49523d0-5543-11ec-3eec-fbe419cf66af
begin

# This cell computes the estimates for column (3) and stores the results.

res3 = estimate(dt,3)
param3 = res3[1]
like3  = res3[2]
conv3  = res3[3]

end;

# ╔═╡ a3f7d260-5543-11ec-3bd8-3d1836a9524e
begin

# This cell computes the estimates for column (4) and stores the results.

res4 = estimate(dt,4)
param4 = res4[1]
like4  = res4[2]
conv4  = res4[3]

end;

# ╔═╡ 7ccc8cbe-5544-11ec-3b68-31c23aaacd87
begin

# unpack the results:

# The params vectors we obtained are a list of lists where the first element is a list containing the estimates for the first individual in the dataset, the second element is a list containing the estimates for the second individual in the dataset, etc.

# For column 1:
b1  = [item[1] for item in param1]; bh1 = [item[2] for item in param1]
d1  = [item[3] for item in param1]; g1  = [item[4] for item in param1]
p1  = [item[5] for item in param1]; s1  = [item[6] for item in param1]

# For column 2:
b2  = [item[1] for item in param2]; bh2 = [item[2] for item in param2]
d2  = [item[3] for item in param2]; g2  = [item[4] for item in param2]
p2  = [item[5] for item in param2]; s2  = [item[6] for item in param2]

# For column 3:
b3  = [item[1] for item in param3]; bh3 = [item[2] for item in param3]
d3  = [item[3] for item in param3]; g3  = [item[4] for item in param3]
p3  = [item[5] for item in param3]; s3  = [item[6] for item in param3]

# For column 4:
b4  = [item[1] for item in param4]; bh4 = [item[2] for item in param4]
d4  = [item[3] for item in param4]; g4  = [item[4] for item in param4]
p4  = [item[5] for item in param4]; s4  = [item[6] for item in param4]
a4 =  [item[7] for item in param4]

end;

# ╔═╡ 8617c2c0-55d2-11ec-007f-315cd7023c9f
# 4. Print and Save Estimation Results

# From the initial 100 subjects, the authors drop those subjects whose estimates their Stata algorithm was unable to find, took longer than 200 iterations to converge, or were outliers (identified by a grubbs test). Since we are using a different programming language and a different minimization algorithm (which, in our case, is gradient-free), we managed to find estimates for all 100 individuals. Below, we show results for both the subset of individuals considered by the authors (which we identified by running their "03MergeIndMLEAndConstructMainSample.do" file) and the whole sample.

# First, we show results for the subset of individuals considered by the authors.

# ╔═╡ 861774a0-55d2-11ec-23d0-2feb0b638bc2
begin

# Load the dataset containing the wid of the individuals to keep in the different specifications. We created this csv file by running their "03MergeIndMLEAndConstructMainSample"

# 4 columns, each containing individuals to keep in each specification
ind_keep = readdlm(raw"../input/ind_to_keep.csv",',',header=true)

# keep parameters only for individuals in ind_keep. First create a dataframe containing the estimates for column 1, 2, 3, 4 and then remove individuals to drop

# For column 1:
mc1 = DataFrame()
mc1.wid=unique(dt.wid); mc1.b=b1; mc1.bh=bh1; mc1.d=d1; mc1.g=g1;
dc1 = mc1[[i for i=1:length(mc1.wid) if mc1.wid[i] in(ind_keep[1][:,1])],:]

# For column 2:
mc2 = DataFrame()
mc2.wid=unique(dt.wid); mc2.b=b2; mc2.bh=bh2; mc2.d=d2; mc2.g=g2;
dc2 = mc2[[i for i=1:length(mc2.wid) if mc2.wid[i] in(ind_keep[1][:,2])],:]

# For column 3:
mc3 = DataFrame()
mc3.wid=unique(dt.wid); mc3.b=b3; mc3.bh=bh3; mc3.d=d3; mc3.g=g3;
dc3 = mc3[[i for i=1:length(mc3.wid) if mc3.wid[i] in(ind_keep[1][:,3])],:]

# For column 4:
mc4 = DataFrame()
mc4.wid=unique(dt.wid); mc4.b=b4; mc4.bh=bh4; mc4.d=d4; mc4.g=g4; mc4.a=a4
dc4 = mc4[[i for i=1:length(mc4.wid) if mc4.wid[i] in(ind_keep[1][:,4])],:]

end;

# ╔═╡ a49dbdd0-55dc-11ec-1e48-21d25efc1b43
begin

# Print the results using the authors' subset of individuals

pnames1 = ["mean(β)", "median(β)", "sd(β)", "mean(β_h)", "median(β_h)", "sd(β_h)",               "mean(δ)", "median(δ)", "sd(δ)", "mean(γ)", "median(γ)", "sd(γ)","mean(α)",           "median(α)", "sd(α)", "P[β]<1", "P[β_h]<1", "r(β,β_h)","p-value r(β,β_h)",
          "Observations"]

at1 = DataFrame()
at1.parameters = pnames1
at1.Primary_Estimation = round.([mean(dc1.b);median(dc1.b);std(dc1.b);mean(dc1.bh);median(dc1.bh);std(dc1.bh);mean(dc1.d);median(dc1.d);std(dc1.d);mean(dc1.g);median(dc1.g);std(dc1.g);0.0;0.0;0.0;mean(dc1.b.<1);mean(dc1.bh.<1);cor(dc1.b,dc1.bh);pvalue(CorrelationTest(dc1.b,dc1.bh));length(dc1.wid)],digits=3)
at1
at1.Early_Decisions = round.([mean(dc2.b);median(dc2.b);std(dc2.b);mean(dc2.bh);median(dc2.bh);std(dc2.bh);mean(dc2.d);median(dc2.d);std(dc2.d);mean(dc2.g);median(dc2.g);std(dc2.g);0.0;0.0;0.0;mean(dc2.b.<1);mean(dc2.bh.<1);cor(dc2.b,dc2.bh);pvalue(CorrelationTest(dc2.b,dc2.bh));length(dc2.wid)],digits=3)
at1
at1.Later_Decisions = round.([mean(dc3.b);median(dc3.b);std(dc3.b);mean(dc3.bh);median(dc3.bh);std(dc3.bh);mean(dc3.d);median(dc3.d);std(dc3.d);mean(dc3.g);median(dc3.g);std(dc3.g);0.0;0.0;0.0;mean(dc3.b.<1);mean(dc3.bh.<1);cor(dc3.b,dc3.bh);pvalue(CorrelationTest(dc3.b,dc3.bh));length(dc3.wid)],digits=3)
at1
at1.Proj_Bias = round.([mean(dc4.b);median(dc4.b);std(dc4.b);mean(dc4.bh);median(dc4.bh);std(dc4.bh);mean(dc4.d);median(dc4.d);std(dc4.d);mean(dc4.g);median(dc4.g);std(dc4.g);mean(dc4.a);median(dc4.a);std(dc4.a);mean(dc4.b.<1);mean(dc4.bh.<1);cor(dc4.b,dc4.bh);pvalue(CorrelationTest(dc4.b,dc4.bh));length(dc4.wid)],digits=3)
at1

end

# ╔═╡ 9087d8e0-55e0-11ec-35df-21166bf459e5
begin

# Second, we present results including our estimates for those subjects the Stata algorithm was unable to converge in less than 200 iterations. As the authors do, we remove outliers using a grubbs test with alpha at 1%. There is no package in Julia that computes the grubbs test for us, so we need to write a function ourselves:

function grubbs_test(y::Vector{Float64}, alpha)

	# compute test statistics
	x = copy(y)
	avg_x = mean(x)
	abs_val_minus_avg = abs.(x .- avg_x)
	max_dev = maximum(abs_val_minus_avg)
	std_x = std(x)
	G = max_dev / std_x
	index = findall(x->x==max_dev, abs_val_minus_avg)

	# compute critical value

	t_dist = quantile.(TDist(length(x)-2), alpha / (2*length(x)))
	numerator = (length(x)-1) * sqrt(t_dist^2)
	denominator = sqrt(length(x)) * sqrt(length(x)-2+t_dist^2)
	critical_value = numerator / denominator

	if G > critical_value deleteat!(x,index) end   # remove if outlier

	return x

end

# We now write a function for the Generalized ESD test for outliers. We basically reiterate the function grubbs_test until no outliers are detected

function esd_test(x::Vector{Float64}, alpha)

	diff = 1

	while diff != 0
		oldlength = length(x)
		x = grubbs_test(x,alpha)
		newlength = length(x)
		diff = oldlength - newlength
	end

	return x

end

end

# ╔═╡ 08059440-55f3-11ec-0d98-438f14370058
begin

# Apply ESD test

beta1 = esd_test(b1, 0.01);  betahat1 = esd_test(bh1, 0.01)
delta1 = esd_test(d1, 0.01); gamma1 = esd_test(g1, 0.01)
beta2 = esd_test(b2, 0.01);  betahat2 = esd_test(bh2, 0.01)
delta2 = esd_test(d2, 0.01); gamma2 = esd_test(g2, 0.01)
beta3 = esd_test(b3, 0.01);  betahat3 = esd_test(bh3, 0.01)
delta3 = esd_test(d3, 0.01); gamma3 = esd_test(g3, 0.01)
beta4 = esd_test(b4, 0.01);  betahat4 = esd_test(bh4, 0.01)
delta4 = esd_test(d4, 0.01); gamma4 = esd_test(g4, 0.01)
alpha4 = esd_test(a4, 0.01)

# we keep only those individuals for which none of the parameters are removed by the grubbs test

# For column 1:

mdc1 = mc1[[i for i=1:length(mc1.wid) if mc1.b[i] in(beta1)],:]
mdc1 = mdc1[[i for i=1:length(mdc1.wid) if mdc1.bh[i] in(betahat1)],:]
mdc1 = mdc1[[i for i=1:length(mdc1.wid) if mdc1.d[i] in(delta1)],:]
mdc1 = mdc1[[i for i=1:length(mdc1.wid) if mdc1.g[i] in(gamma1)],:]

# For column 2:

mdc2 = mc2[[i for i=1:length(mc2.wid) if mc2.b[i] in(beta2)],:]
mdc2 = mdc2[[i for i=1:length(mdc2.wid) if mdc2.bh[i] in(betahat2)],:]
mdc2 = mdc2[[i for i=1:length(mdc2.wid) if mdc2.d[i] in(delta2)],:]
mdc2 = mdc2[[i for i=1:length(mdc2.wid) if mdc2.g[i] in(gamma2)],:]

# For column 3:

mdc3 = mc3[[i for i=1:length(mc3.wid) if mc3.b[i] in(beta3)],:]
mdc3 = mdc3[[i for i=1:length(mdc3.wid) if mdc3.bh[i] in(betahat3)],:]
mdc3 = mdc3[[i for i=1:length(mdc3.wid) if mdc3.d[i] in(delta3)],:]
mdc3 = mdc3[[i for i=1:length(mdc3.wid) if mdc3.g[i] in(gamma3)],:]

# For column 4:

mdc4 = mc4[[i for i=1:length(mc4.wid) if mc4.b[i] in(beta4)],:]
mdc4 = mdc4[[i for i=1:length(mdc4.wid) if mdc4.bh[i] in(betahat4)],:]
mdc4 = mdc4[[i for i=1:length(mdc4.wid) if mdc4.d[i] in(delta4)],:]
mdc4 = mdc4[[i for i=1:length(mdc4.wid) if mdc4.g[i] in(gamma4)],:]
mdc4 = mdc4[[i for i=1:length(mdc4.wid) if mdc4.a[i] in(alpha4)],:]

end;

# ╔═╡ c4bcd1a0-55f5-11ec-357b-87eef6d8494e
begin

# Print our results

pnames2 = ["mean(β)", "median(β)", "sd(β)", "mean(β_h)", "median(β_h)", "sd(β_h)",               "mean(δ)", "median(δ)", "sd(δ)", "mean(γ)", "median(γ)", "sd(γ)","mean(α)",           "median(α)", "sd(α)", "P[β]<1", "P[β_h]<1", "r(β,β_h)","p-value r(β,β_h)",
          "Observations"]

at2 = DataFrame()
at2.parameters = pnames2
at2.Primary_Estimation = round.([mean(mdc1.b);median(mdc1.b);std(mdc1.b);mean(mdc1.bh);median(mdc1.bh);std(mdc1.bh);mean(mdc1.d);median(mdc1.d);std(mdc1.d);mean(mdc1.g);median(mdc1.g);std(mdc1.g);0.0;0.0;0.0;mean(mdc1.b.<1);mean(mdc1.bh.<1);cor(mdc1.b,mdc1.bh);pvalue(CorrelationTest(mdc1.b,mdc1.bh));length(mdc1.wid)],digits=3)
at2.Early_Decisions = round.([mean(mdc2.b);median(mdc2.b);std(mdc2.b);mean(mdc2.bh);median(mdc2.bh);std(mdc2.bh);mean(mdc2.d);median(mdc2.d);std(mdc2.d);mean(mdc2.g);median(mdc2.g);std(mdc2.g);0.0;0.0;0.0;mean(mdc2.b.<1);mean(mdc2.bh.<1);cor(mdc2.b,mdc2.bh);pvalue(CorrelationTest(mdc2.b,mdc2.bh));length(mdc2.wid)],digits=3)
at2.Later_Decisions = round.([mean(mdc3.b);median(mdc3.b);std(mdc3.b);mean(mdc3.bh);median(mdc3.bh);std(mdc3.bh);mean(mdc3.d);median(mdc3.d);std(mdc3.d);mean(mdc3.g);median(mdc3.g);std(mdc3.g);0.0;0.0;0.0;mean(mdc3.b.<1);mean(mdc3.bh.<1);cor(mdc3.b,mdc3.bh);pvalue(CorrelationTest(mdc3.b,mdc3.bh));length(mdc3.wid)],digits=3)
at2.Proj_Bias = round.([mean(mdc4.b);median(mdc4.b);std(mdc4.b);mean(mdc4.bh);median(mdc4.bh);std(mdc4.bh);mean(mdc4.d);median(mdc4.d);std(mdc4.d);mean(mdc4.g);median(mdc4.g);std(mdc4.g);mean(mdc4.a);median(mdc4.a);std(mdc4.a);mean(mdc4.b.<1);mean(mdc4.bh.<1);cor(mdc4.b,mdc4.bh);pvalue(CorrelationTest(mdc4.b,mdc4.bh));length(mdc4.wid)],digits=3)
at2

end

# ╔═╡ c072cbde-55f5-11ec-1a99-c17e1667c4e2
begin

# Save our raw results

table2_julia = DataFrame()
table2_julia.beta1 = b1
table2_julia.betahat1 = bh1
table2_julia.delta1 = d1
table2_julia.gamma1 = g1
table2_julia.phi1 = p1
table2_julia.sigma1 = s1
table2_julia.beta2 = b2
table2_julia.betahat2 = bh2
table2_julia.delta2 = d2
table2_julia.gamma2 = g2
table2_julia.phi2 = p2
table2_julia.sigma2 = s2
table2_julia.beta3 = b3
table2_julia.betahat3 = bh3
table2_julia.delta3 = d3
table2_julia.gamma3 = g3
table2_julia.phi3 = p3
table2_julia.sigma3 = s3
table2_julia.beta4 = b4
table2_julia.betahat4 = bh4
table2_julia.delta4 = d4
table2_julia.gamma4 = g4
table2_julia.phi4 = p4
table2_julia.sigma4 = s4
table2_julia.alpha4 = a4

CSV.write(raw"../output/table2_julia.csv", table2_julia)

end

# ╔═╡ fe89b822-268f-11ec-1759-67648ab10a25

# ╔═╡ fe75baf0-268f-11ec-3b55-15fbd2599c0f

# ╔═╡ 4c4dd50e-2757-11ec-3b42-f76583bd9178

# ╔═╡ e5515240-2758-11ec-0aa4-c1dba4b2b512

# ╔═╡ 62c2f7fe-2759-11ec-0cd8-33f81746ff23

# ╔═╡ ca7ac9a2-2759-11ec-1411-afe4da3e4d7d

# ╔═╡ Cell order:
# ╠═d08aaa10-268f-11ec-067f-4322e82ba779
# ╠═fee9d980-268f-11ec-3901-3705290b7c24
# ╠═f8c1b410-51b4-11ec-1d65-11df2949bf23
# ╠═0de27c80-51b5-11ec-05e9-a5b37469b0c7
# ╠═feb1b280-268f-11ec-1574-e5037552b641
# ╠═fe9ddc60-268f-11ec-112a-59ea860b2c81
# ╠═74270f10-51b5-11ec-2edd-af650c9bec6b
# ╠═ce679cb0-51b5-11ec-39ba-3f65d7b1ce33
# ╠═a5d0925e-5535-11ec-06ca-85bfd5b5081b
# ╠═f2177020-5536-11ec-0db6-1b2a59aa1c0e
# ╠═99348ad0-5543-11ec-274d-1d59ea3e02b0
# ╠═a49523d0-5543-11ec-3eec-fbe419cf66af
# ╠═a3f7d260-5543-11ec-3bd8-3d1836a9524e
# ╠═7ccc8cbe-5544-11ec-3b68-31c23aaacd87
# ╠═8617c2c0-55d2-11ec-007f-315cd7023c9f
# ╠═861774a0-55d2-11ec-23d0-2feb0b638bc2
# ╠═a49dbdd0-55dc-11ec-1e48-21d25efc1b43
# ╠═9087d8e0-55e0-11ec-35df-21166bf459e5
# ╠═08059440-55f3-11ec-0d98-438f14370058
# ╠═c4bcd1a0-55f5-11ec-357b-87eef6d8494e
# ╠═c072cbde-55f5-11ec-1a99-c17e1667c4e2
# ╟─fe89b822-268f-11ec-1759-67648ab10a25
# ╟─fe75baf0-268f-11ec-3b55-15fbd2599c0f
# ╟─4c4dd50e-2757-11ec-3b42-f76583bd9178
# ╟─e5515240-2758-11ec-0aa4-c1dba4b2b512
# ╟─62c2f7fe-2759-11ec-0cd8-33f81746ff23
# ╟─ca7ac9a2-2759-11ec-1411-afe4da3e4d7d


text:
module TestSCS

using Test
using MathOptInterface
import SCS

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

test_DirectSolver() = _test_runtests(SCS.DirectSolver)

test_IndirectSolver() = _test_runtests(SCS.IndirectSolver)

function _test_runtests(linear_solver)
    optimizer = SCS.Optimizer()
    MOI.set(
        optimizer,
        MOI.RawOptimizerAttribute("linear_solver"),
        linear_solver,
    )
    MOI.set(optimizer, MOI.RawOptimizerAttribute("eps"), 1e-6)
    MOI.set(optimizer, MOI.Silent(), true)
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            optimizer,
        ),
        Float64,
    )
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            atol = 1e-4,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
                MOI.ConstraintName,
                MOI.VariableName,
                MOI.ObjectiveBound,
            ],
        ),
        exclude = String[
            # Expected test failures:
            #   ArgumentError: The number of constraints must be greater than 0
            "test_attribute_RawStatusString",
            "test_attribute_SolveTimeSec",
            "test_objective_ObjectiveFunction_blank",
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
            #   Problem is a nonconvex QP
            "test_basic_ScalarQuadraticFunction_EqualTo",
            "test_basic_ScalarQuadraticFunction_GreaterThan",
            "test_basic_ScalarQuadraticFunction_Interval",
            "test_basic_VectorQuadraticFunction_",
            "test_quadratic_SecondOrderCone_basic",
            "test_quadratic_nonconvex_",
            #   power cone error, values must be in [-1,1]
            "test_conic_DualPowerCone_VectorOfVariables",
            "test_conic_DualPowerCone_VectorAffineFunction",
            "test_conic_PowerCone_VectorAffineFunction",
            "test_conic_PowerCone_VectorOfVariables",
            #   MathOptInterface.jl issue #1431
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
        ],
    )
    return
end

function test_RawOptimizerAttribute()
    model = SCS.Optimizer()
    MOI.set(model, MOI.RawOptimizerAttribute("eps"), 1.0)
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 1.0
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 1.0
    MOI.set(model, MOI.RawOptimizerAttribute("eps"), 2.0)
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 2.0
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps")) == 2.0
end

function test_constrained_variables()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            SCS.Optimizer(),
        ),
        Float64,
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.PositiveSemidefiniteConeTriangle,
    )
    x = MOI.add_variables(model, 6)
    f = MOI.VectorOfVariables(x)
    s = MOI.PositiveSemidefiniteConeTriangle(3)
    @test isa(
        MOI.add_constraint(model, f, s),
        MOI.ConstraintIndex{typeof(f),typeof(s)},
    )
    return
end

end  # module

TestSCS.runtests()

ext: jl


text:
function getgensym!(defined::Dict{Symbol,Symbol}, s::Symbol)
  snew = get!(defined, s) do
    gensym(s)
  end
  snew ≢ s && (defined[snew] = snew)
  return snew
end

extractargs!(arguments::Vector{Symbol}, defined::Dict{Symbol,Symbol}, sym, mod) = nothing
# function extractargs!(arguments::Vector{Symbol}, defined::Dict{Symbol,Symbol}, sym::Symbol, mod)
#   if ((sym ∉ keys(defined)) && sym ∉ (:nothing, :(+), :(*), :(-), :(/), :(÷), :(<<), :(>>), :(>>>), :zero, :one)) && !Base.isdefined(mod, sym)
#     @show getgensym!(defined, sym)
#     @assert false
#     push!(arguments, getgensym!(defined, sym))
#   end
#   nothing
# end
function define_tup!(arguments::Vector{Symbol}, defined::Dict{Symbol,Symbol}, ex::Expr, mod)
  for (i, a) ∈ enumerate(ex.args)
    if a isa Symbol
      ex.args[i] = getgensym!(defined, a)
    elseif Meta.isexpr(a, :tuple)
      define_tup!(defined, a)
    elseif Meta.isexpr(a, :ref)
      extractargs!(arguments, defined, a, mod)
    else
      throw("Don't know how to handle:\n $a")
    end
  end
end
function define1!(arguments::Vector{Symbol}, defined::Dict{Symbol,Symbol}, x::Vector{Any}, mod)
  s = x[1]
  if s isa Symbol
    x[1] = getgensym!(defined, s)
  else
    define_tup!(arguments, defined, s::Expr, mod)
  end
end
function define_induction_variables!(arguments::Vector{Symbol}, defined::Dict{Symbol,Symbol}, ex::Expr, mod) # add `i` in `for i ∈ looprange` to `defined`
  ex.head === :for || return
  loops = ex.args[1]
  if loops.head === :block
    for loop ∈ loops.args
      define1!(arguments, defined, loop.args, mod)
    end
  else
    define1!(arguments, defined, loops.args, mod)
  end
end

function extractargs_equal!(arguments::Vector{Symbol}, defined::Dict{Symbol,Symbol}, args::Vector{Any}, mod)
  arg1 = first(args)
  if arg1 isa Symbol
    args[1] = getgensym!(defined, arg1)
  elseif Meta.isexpr(arg1, :tuple)
    define_tup!(arguments, defined, arg1, mod)
  end
  nothing
end
function must_add_sym(defined::Dict{Symbol,Symbol}, arg::Symbol, mod)
  ((arg ∉ keys(defined)) && arg ∉ (:nothing, :(+), :(*), :(-), :(/), :(÷), :(<<), :(>>), :(>>>), :zero, :one)) && !Base.isdefined(mod, arg)
end
function get_sym!(defined::Dict{Symbol,Symbol}, arguments::Vector{Symbol}, arg::Symbol, mod)
  if must_add_sym(defined, arg, mod)
    # @show getgensym!(defined, sym)
    # @assert false
    push!(arguments, arg)
    getgensym!(defined, arg)
  else
    get(defined, arg, arg)
  end
end
function extractargs!(arguments::Vector{Symbol}, defined::Dict{Symbol,Symbol}, expr::Expr, mod)
  define_induction_variables!(arguments, defined, expr, mod)
  head = expr.head
  args = expr.args
  startind = 1
  if head === :call
    startind = 2
  elseif head === :(=)
    extractargs_equal!(arguments, defined, args, mod)
  elseif head ∈ (:inbounds, :loopinfo)#, :(->))
    return
  elseif head === :(.)
    arg1 = args[1]
    if arg1 isa Symbol
      args[1] = get_sym!(defined, arguments, arg1, mod)
    else
      extractargs!(arguments, defined, arg1, mod)
    end
  elseif head === :(->)
    td = copy(defined)
    define1!(arguments::Vector{Symbol}, td, args, mod)
    extractargs!(arguments, td, args[1], mod)
    extractargs!(arguments, td, args[2], mod)
    return
  elseif (head === :local) || (head === :global)
    for (i,arg) in enumerate(args)
      if Meta.isexpr(arg, :(=))
        extractargs_equal!(arguments, defined, arg.args, mod)
        args = arg.args
        startind = 2
      else
        args[i] = getgensym!(defined, arg)
        return
      end
    end
  elseif head === :kw
    return
  end
  for i ∈ startind:length(args)
    argᵢ = args[i]
    (head === :ref && ((argᵢ === :end) || (argᵢ === :begin))) && continue
    if argᵢ isa Symbol
      args[i] = get_sym!(defined, arguments, argᵢ, mod)
    elseif argᵢ isa Expr
      extractargs!(arguments, defined, argᵢ, mod)
    else
      extractargs!(arguments, defined, argᵢ, mod)
    end
  end
  return
end

function symbolsubs(e::Expr, old::Symbol, new::Symbol)
  return Expr(e.head, (symbolsubs(a, old, new) for a in e.args)...)
end
function symbolsubs(e::Symbol, old::Symbol, new::Symbol)
  e==old ? new : e
end
symbolsubs(e, old::Symbol, new::Symbol) = e

struct NoLoop end
Base.iterate(::NoLoop) = (NoLoop(), NoLoop())
Base.iterate(::NoLoop, ::NoLoop) = nothing
@inline splitloop(x) = NoLoop(), x, CombineIndices()
struct CombineIndices end
@inline splitloop(x::AbstractRange) = NoLoop(), x, CombineIndices()
@inline function splitloop(x::CartesianIndices)
  axes = x.indices
  CartesianIndices(Base.front(axes)), last(axes), CombineIndices()
end
@inline function splitloop(x::AbstractArray)
  inds = eachindex(x)
  inner, outer = splitloop(inds)
  inner, outer, x
end
struct TupleIndices end
@inline function splitloop(x::Base.Iterators.ProductIterator{Tuple{T1,T2}}) where {T1,T2}
  iters = x.iterators
  iters[1], iters[2], TupleIndices()
end
@inline function splitloop(x::Base.Iterators.ProductIterator{<:Tuple{Vararg{Any,N}}}) where {N}
  iters = x.iterators
  Base.front(iters), iters[N], TupleIndices()
end
combine(::CombineIndices, ::NoLoop, x) = x
combine(::CombineIndices, I::CartesianIndex, j) = CartesianIndex((I.I..., j))
combine(::TupleIndices, i::Tuple, j) = (i..., j)
combine(::TupleIndices, i::Number, j) = (i, j)

Base.@propagate_inbounds combine(x::AbstractArray, I, j) = x[combine(CombineIndices(), I, j)]
Base.@propagate_inbounds combine(x::AbstractArray, ::NoLoop, j) = x[j]

static_literals!(s::Symbol) = s
function static_literals!(q::Expr)
  for (i,ex) ∈ enumerate(q.args)
    if ex isa Integer
      q.args[i] = StaticInt(ex)
    elseif ex isa Expr
      static_literals!(ex)
    end
  end
  q
end
function maybestatic!(_expr)::Expr
  _expr isa Expr || return esc(_expr)
  expr::Expr = _expr
  if expr.head === :call
    f = first(expr.args)
    if f === :length
      expr.args[1] = static_length
    elseif f === :size && length(expr.args) == 3
      i = expr.args[3]
      if i isa Integer
        expr.args[1] = size
        expr.args[3] = StaticInt(i)
      end
    else
      static_literals!(expr)
    end
  end
  esc(expr)
end
function enclose(exorig::Expr, reserve_per, minbatchsize, per::Symbol, threadlocal, mod)
  Meta.isexpr(exorig, :for, 2) || throw(ArgumentError("Expression invalid; should be a for loop."))
  ex = copy(exorig)
  loop_sym = Symbol("##LOOP##")
  loopstart = Symbol("##LOOPSTART##")
  loop_step = Symbol("##LOOP_STEP##")
  loop_stop = Symbol("##LOOP_STOP##")
  iter_leng = Symbol("##ITER_LENG##")
  loop_offs = Symbol("##LOOPOFFSET##")
  innerloop = Symbol("##inner##loop##")
  rcombiner = Symbol("##split##recombined##")
  threadlocal_var = Symbol("threadlocal")

  # arguments = Symbol[]#loop_offs, loop_step]
  arguments = Symbol[innerloop, rcombiner]#loop_offs, loop_step]
  defined = Dict{Symbol,Symbol}(loop_offs => loop_offs, loop_step => loop_step)
  threadlocal_var_gen = getgensym!(defined, threadlocal_var)
  define_induction_variables!(arguments, defined, ex, mod)
  firstloop = ex.args[1]
  if firstloop.head === :block
    secondaryloopsargs = firstloop.args[2:end]
    firstloop = firstloop.args[1]
  else
    secondaryloopsargs = Any[]
  end
  loop = firstloop.args[2]
  # @show ex loop
  firstlooprange = Expr(:call, GlobalRef(Base, :(:)), loopstart, loop_step, loop_stop)
  body = ex.args[2]
  if length(secondaryloopsargs) == 1
    body = Expr(:for, only(secondaryloopsargs), body)
  elseif length(secondaryloopsargs) > 1
    sl = Expr(:block); append!(sl.args, secondaryloopsargs)
    body = Expr(:for, sl, body)
  end
  fla1 = firstloop.args[1]
  excomb = if fla1 isa Symbol
    fla1 = getgensym!(defined, fla1)
    quote
      # for $(firstloop.args[1]) in
      for var"##outer##" in $firstlooprange, var"##inner##" in $innerloop
        $fla1 = $combine($rcombiner, var"##inner##", var"##outer##")
        $body
      end
    end
  else
    @assert fla1 isa Expr
    for i in eachindex(fla1.args)
      fla1.args[i] = getgensym!(defined, fla1.args[i])
    end
    quote
      # for $(firstloop.args[1]) in
      for var"##outer##" in $firstlooprange, var"##inner##" in $innerloop
        $fla1 = $combine($rcombiner, var"##inner##", var"##outer##")
        $body
      end
    end
  end
  if ex.args[1].head === :block
    for i ∈ 2:length(ex.args[1].args)
      extractargs!(arguments, defined, ex.args[1].args[i], mod)
    end
  end
  for i ∈ 2:length(ex.args)
    extractargs!(arguments, defined, ex.args[i], mod)
  end
  # @show ex.args[1] firstloop body
  # if length(ex.args[
  # ex = quote
  #   # for $(firstloop.args[1]) in
  #   for var"##outer##" in $firstlooprange, var"##inner##" in $innerloop
  #     $(firstloop.args[1]) = $combine($rcombiner, var"##inner##", var"##outer##")
  #     $body
  #   end
  # end
  # typexpr_incomplete is missing `funcs`
  # outerloop = Symbol("##outer##")
  q = quote
    $(esc(innerloop)), $loop_sym, $(esc(rcombiner)) = $splitloop($(maybestatic!(loop)))
    # $loop_sym = $(maybestatic!(loop))
    $iter_leng = $static_length($loop_sym)
    $loop_step = $static_step($loop_sym)
    $loop_offs = $static_first($loop_sym)
  end
  threadtup = Expr(:tuple, iter_leng)
  num_thread_expr = Expr(:call, num_threads)
  if per === :core
    num_thread_expr = Expr(:call, min, num_thread_expr, Expr(:call, num_cores))
  end
  if minbatchsize isa Integer && minbatchsize ≤ 1
    # if reserve_per ≤ 0
      push!(threadtup.args, :(min($iter_leng, $num_thread_expr)))
    # else
    #   push!(threadtup.args, :(min($iter_leng, cld($num_thread_expr, $reserve_per))), reserve_per)
    # end
  else
    il = :(div($iter_leng, $(minbatchsize isa Int ? StaticInt(minbatchsize) : esc(minbatchsize))))
    # if reserve_per ≤ 0
      push!(threadtup.args, :(min($il, $num_thread_expr)))
    # else
    #   push!(threadtup.args, :(min($il, cld($num_thread_expr, $reserve_per))), reserve_per)
    # end
  end
  closure = Symbol("##closure##")
  threadlocal, threadlocal_type = threadlocal
  threadlocal_var_single = gensym(threadlocal_var)
  q_single = symbolsubs(exorig, threadlocal_var, threadlocal_var_single)
  threadlocal_init_single   = threadlocal == :() ? :() : :( $threadlocal_var_single = $threadlocal )
  threadlocal_repack_single = threadlocal == :() ? :() : :( $threadlocal_var_single )
  threadlocal_init1         = threadlocal == :() ? :() : :( $threadlocal_var = Vector{$threadlocal_type}(undef, 0) )
  threadlocal_init2         = threadlocal == :() ? :() : :( resize!($(esc(threadlocal_var)),max(1,$(threadtup.args[2]))) )
  threadlocal_get           = threadlocal == :() ? :() : :( $threadlocal_var_gen = $threadlocal::$threadlocal_type )
  threadlocal_set           = threadlocal == :() ? :() : :( $threadlocal_var[var"##THREAD##"] = $threadlocal_var_gen )
  push!(q.args, threadlocal_init2)
  args = Expr(:tuple, Symbol("##LOOPOFFSET##"), Symbol("##LOOP_STEP##"))
  closure_args = if threadlocal == :()
    :($args, var"##SUBSTART##"::Int, var"##SUBSTOP##"::Int)
  else
    :($args, var"##SUBSTART##"::Int, var"##SUBSTOP##"::Int, var"##THREAD##"::Int)
  end
  closureq = quote
    $closure = let
      @inline $closure_args -> begin
        var"##LOOPSTART##" = var"##SUBSTART##" * var"##LOOP_STEP##" + var"##LOOPOFFSET##" - var"##LOOP_STEP##"
        var"##LOOP_STOP##" = var"##SUBSTOP##" * var"##LOOP_STEP##" + var"##LOOPOFFSET##" - var"##LOOP_STEP##"
        $threadlocal_get
        @inbounds begin
          $excomb
        end
        $threadlocal_set
        nothing
      end
    end
  end
  push!(q.args, esc(closureq))
  batchcall = Expr(:call, batch, esc(closure), threadtup, Symbol("##LOOPOFFSET##"), Symbol("##LOOP_STEP##"))
  for a ∈ arguments
    push!(args.args, get(defined,a,a))
    push!(batchcall.args, esc(a))
  end
  if threadlocal != :()
    push!(batchcall.args, Expr(:kw, :threadlocal, true))
  end
  push!(q.args, batchcall)
  quote
    if $num_threads() == 1
      single_thread_result = begin
        $(esc(threadlocal_init_single)) # Initialize threadlocal storage
        $(esc(q_single))
        $(esc(threadlocal_repack_single))
      end
      $(esc(threadlocal_var)) = [single_thread_result] # Put the single-thread threadlocal storage in a single-element Vector
    else
      $(esc(threadlocal_init1))
      let
        $q
      end
    end
  end
end

"""
    @batch for i in Iter; ...; end

Evaluate the loop on multiple threads.

    @batch minbatch=N for i in Iter; ...; end

Create a thread-local storage used in the loop.

    @batch localthread=init() for i in Iter; ...; end

The `init` function will be called at the start at each thread. `localthread` will
refer to storage local for the thread. At the end of the loop, a `localthread`
vector containing all the thread-local values will be available. A type can be specified
with `threadlocal=init()::Type`.

Evaluate at least N iterations per thread. Will use at most `length(Iter) ÷ N` threads.

    @batch per=core for i in Iter; ...; end
    @batch per=thread for i in Iter; ...; end

Use at most 1 thread per physical core, or 1 thread per CPU thread, respectively.
One thread per core will mean less threads competing for the cache, while (for example) if
there are two hardware threads per physical core, then using each thread means that there
are two independent instruction streams feeding the CPU's execution units. When one of
these streams isn't enough to make the most of out of order execution, this could increase
total throughput.

Which performs better will depend on the workload, so if you're not sure it may be worth
benchmarking both.

LoopVectorization.jl currently only uses up to 1 thread per physical core. Because there
is some overhead to switching the number of threads used, `per=core` is `@batch`'s default,
so that `Polyester.@batch` and `LoopVectorization.@tturbo` work well together by default.

You can pass both `per=(core/thread)` and `minbatch=N` options at the same time, e.g.

    @batch per=thread minbatch=2000 for i in Iter; ...; end
    @batch minbatch=5000 per=core   for i in Iter; ...; end
"""
macro batch(ex)
  enclose(macroexpand(__module__, ex), 0, 1, :core, (:(),:Any), __module__)
end
function interpret_kwarg(arg, reserve_per = 0, minbatch = 1, per = :core, threadlocal = (:(),:Any))
  a = arg.args[1]
  v = arg.args[2]
  if a === :reserve
    @assert v ≥ 0
    reserve_per = v
  elseif a === :minbatch
    minbatch = v
  elseif a === :per
    per = v::Symbol
    @assert (per === :core) | (per === :thread)
  elseif a === :threadlocal
    if isa(v,Expr) && v.head==:(::)
      threadlocal = v.args
    else
      threadlocal = (v, :Any)
    end
  else
    throw(ArgumentError("kwarg $(a) not recognized."))
  end
  reserve_per, minbatch, per, threadlocal
end
macro batch(arg1, ex)
  reserve, minbatch, per, threadlocal = interpret_kwarg(arg1)
  enclose(macroexpand(__module__, ex), reserve, minbatch, per, threadlocal, __module__)
end
macro batch(arg1, arg2, ex)
  reserve, minbatch, per, threadlocal = interpret_kwarg(arg1)
  reserve, minbatch, per, threadlocal = interpret_kwarg(arg2, reserve, minbatch, per, threadlocal)
  enclose(macroexpand(__module__, ex), reserve, minbatch, per, threadlocal, __module__)
end
macro batch(arg1, arg2, arg3, ex)
  reserve, minbatch, per, threadlocal = interpret_kwarg(arg1)
  reserve, minbatch, per, threadlocal = interpret_kwarg(arg2, reserve, minbatch, per, threadlocal)
  reserve, minbatch, per, threadlocal = interpret_kwarg(arg3, reserve, minbatch, per, threadlocal)
  enclose(macroexpand(__module__, ex), reserve, minbatch, per, threadlocal, __module__)
end


text:
using Statistics
function mean(ds::Array{Dipole})

    mean_x = Statistics.mean([d.x for d in ds])
    mean_y = Statistics.mean([d.y for d in ds])
    mean_z = Statistics.mean([d.z for d in ds])
    mean_s = Statistics.mean([d.size for d in ds])

    Dipole(ds[1].coord_system, mean_x, mean_y, mean_z, 0, 0, 0, 0, 0, mean_s)
end

function std(ds::Array{Dipole})

    std_x = Statistics.std([d.x for d in ds])
    std_y = Statistics.std([d.y for d in ds])
    std_z = Statistics.std([d.z for d in ds])
    std_s = Statistics.std([d.size for d in ds])

    Dipole(ds[1].coord_system, std_x, std_y, std_z, 0, 0, 0, 0, 0, std_s)
end


text:
################################################################################

@testset "Process Hidden Markov Model" begin

  @test process()
end

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

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
using Documenter, DocumenterCitations, DeconvOptim

cite_bib = CitationBibliography(joinpath(@__DIR__, "../paper/ref.bib"))

DocMeta.setdocmeta!(DeconvOptim, :DocTestSetup, :(using DeconvOptim); recursive=true)
makedocs(cite_bib, modules=[DeconvOptim],
         sitename="DeconvOptim.jl",
         doctest = false,
         pages = Any[
                "DeconvOptim.jl" => "index.md",
                "Workflow" => Any[
                        "workflow/basic_workflow.md",
                        "workflow/changing_regularizers.md",
                        "workflow/changing_loss.md",
                        "workflow/3D_dataset.md",
                        "workflow/cuda.md",
                        "workflow/flexible_invert.md",
                        "workflow/performance_tips.md",
                    ],
                "Background" => Any[
                    "background/physical_background.md",
                    "background/mathematical_optimization.md",
                    "background/loss_functions.md",
                    "background/regularizer.md",
                    ],
                "Function references" => Any[
                    "function_references/deconvolution.md",
                    "function_references/loss.md",
                    "function_references/mapping.md",
                    "function_references/regularizer.md",
                    "function_references/utils.md",
                    "function_references/analysis.md",
                    ],
                "References" => "references.md"
                ],
        )

deploydocs(repo = "github.com/roflmaostc/DeconvOptim.jl.git")


text:
module Lamp

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Lamp" => Ahorn.EntityPlacement(
        Maple.Lamp
    ),
    "Lamp (Broken)" => Ahorn.EntityPlacement(
        Maple.BrokenLamp
    ),
)

function Ahorn.selection(entity::Maple.Lamp)
    sprite = Ahorn.getSprite("scenery/lamp", "Gameplay")
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - floor(Int, sprite.width / 4), y - sprite.height, floor(Int, sprite.width / 2), sprite.height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Lamp, room::Maple.Room)
    sprite = Ahorn.getSprite("scenery/lamp", "Gameplay")

    width = floor(Int, sprite.width / 2)
    height = sprite.height

    broken = get(entity.data, "broken", false)
    Ahorn.drawImage(ctx, sprite, -floor(Int, width / 2), -height, width * broken, 0, width, height)
end

end


text:
include("../common/layers.jl")

using Random

mutable struct SimpleCBOW
    in_layer0
    in_layer1
    out_layer
    loss_layer
    params
    grads
    word_vecs

    function SimpleCBOW(vocab_size, hidden_size)
        self = new()
        V, H = vocab_size, hidden_size

        # 重みの初期化
        W_in = 0.01 .* Random.randn(V, H) # Float32の指定は特にしない
        W_out = 0.01 .* Random.randn(H, V) # Float32の指定は特にしない

        # レイヤの生成
        self.in_layer0 = MatMul(W_in)
        self.in_layer1 = MatMul(W_in)
        self.out_layer = MatMul(W_out)
        self.loss_layer = SoftmaxWithLoss()

        # 全ての重みと勾配を配列にまとめる
        layers = [self.in_layer0, self.in_layer1, self.out_layer, self.loss_layer]
        self.params, self.grads = [], []
        for layer in layers
            push!(self.params, layer.params)
            push!(self.grads, layer.grads)
        end

        # 構造体の変数に単語の分散表現を設定
        self.word_vecs = W_in
        return self
    end
end

function forward!(self::SimpleCBOW, contexts, target)
    # contextsの形状は3次元決め打ち
    h0 = forward!(self.in_layer0, contexts[:, 1, :]) # MatMulのforward!
    h1 = forward!(self.in_layer1, contexts[:, 2, :]) # MatMulのforward!
    h = (h0 .+ h1) .* 0.5
    score = forward!(self.out_layer, h) # MatMulのforward!
    loss = forward!(self.loss_layer, score, target) # SoftmaxWithLossのforward!
    return loss
end

function backward!(self::SimpleCBOW; dout=1)
    ds = backward!(self.loss_layer; dout=dout) # SoftmaxWithLossのbackward!
    da = backward!(self.out_layer, ds) # MatMulのbackward!
    da .*= 0.5
    backward!(self.in_layer1, da) # MatMulのbackward!
    backward!(self.in_layer0, da) # MatMulのbackward!
end


text:
"""
    TFIModel{LT<:AbstractLattice} <: AbstractModel{LT}

Transverse field Ising Model.
"""
struct TFIModel{LT<:AbstractLattice} <: AbstractModel{LT}
    J0::Float64
    Γ::Float64
    temp::Float64
    lattice::LT

    function TFIModel(J0::Real, Γ::Real, temp::Real, lattice::LT) where LT<:AbstractLattice
        new{LT}(Float64(J0), Float64(Γ), Float64(temp), lattice)
    end
end


text:
# Display methods (show & print & writemime) tests
using Base64

@testset "Output" begin

@testset "text/plain" begin
    @test occursin("1.00", stringmime("text/plain", 1USD))
    @test occursin("1", stringmime("text/plain", 1JPY))
    @test occursin("CAD", stringmime("text/plain", Basket([1USD, 1CAD])))
    @test occursin("USD", stringmime("text/plain", Basket([1USD, 1CAD])))
    @test occursin(
        "200.00 EUR",
        stringmime("text/plain", Basket([100USD, 200EUR])))
    @test !occursin(
        "CAD", stringmime("text/plain", Basket([1USD, 1CAD, -1CAD])))
    @test stringmime("text/plain", +USD) == "1.00 USD"
    @test stringmime("text/plain", -USD) == "−1.00 USD"

    @test stringmime("text/plain", majorunit(Monetary{:USD, BigInt, 5})) ==
        "1.00000 USD"
    @test stringmime("text/plain", -Monetary(:JPY, precision=2)) ==
        "−1.00 JPY"
    @test stringmime("text/plain", 7.2512Monetary(:XAU; precision=4)) ==
        "7.2512 XAU"

    # with compact IOContext
    buf = IOBuffer()
    show(IOContext(buf, :compact => true), "text/plain", 1USD)
    @test String(take!(buf)) == "1.0USD"
end

@testset "text/latex" begin
    @test stringmime("text/latex", 100USD) == "\$100.00\\,\\mathrm{USD}\$"
    @test stringmime("text/latex", -100JPY) == "\$-100\\,\\mathrm{JPY}\$"
    @test stringmime("text/latex", 0GBP) == "\$0.00\\,\\mathrm{GBP}\$"
    @test stringmime("text/latex", zero(Monetary{:EUR, Int, 0})) ==
        "\$0\\,\\mathrm{EUR}\$"
end

@testset "text/markdown" begin
    basketstr = stringmime("text/markdown", Basket([USD, 20CAD, -10JPY]))

    @test occursin("`Basket`", basketstr)
    @test occursin("\$3\$-currency", basketstr)
    @test occursin(" - \$-10\\,\\mathrm{JPY}\$", basketstr)
end

@testset "print & show" begin
    @test string(USD) == "1.0USD"
    @test string(0.01USD) == "0.01USD"
    @test string(20JPY) == "20.0JPY"

    # this test is a bit complicated because order is undefined
    basketstr = string(Basket([USD, 20CAD, -10JPY]))
    @test occursin("Basket([", basketstr)
    @test occursin("-10.0JPY", basketstr)
    @test occursin("20.0CAD", basketstr)

    # test compatibility between show & print
    @test sprint(show, 0.02USD) == string(0.02USD)
end

end  # testset output

@testset "format" begin
    # internals
    @test Currencies.CurrencyFormatting.digitseparate(20160408, "'", (3, 3)) == "20'160'408"
    @test Currencies.CurrencyFormatting.digitseparate(12345678, " ", (3, 3)) == "12 345 678"
    @test Currencies.CurrencyFormatting.digitseparate(1, ",", (3, 3)) == "1"
    @test Currencies.CurrencyFormatting.digitseparate(12345678, ",", (3, 2)) == "1,23,45,678"

    # default (finance) style
    @test format(100USD) == "100.00 USD"
    @test format(8050.20USD) == "8050.20 USD"
    @test format(-100USD) == "(100.00) USD"
    @test format(-8050.20EUR) == "(8050.20) EUR"
    @test format(-5JPY) == "(5) JPY"
    @test format(0CAD) == "— CAD"

    # us style
    @test format(1000USD, styles=[:us]) == "USD 1,000.00"
    @test format(1000JPY, styles=[:us]) == "JPY 1,000"
    @test format(19970716.14AUD, styles=[:us]) == "AUD 19,970,716.14"
    @test format(-1100.55USD, styles=[:us]) == "USD −1,100.55"

    # european style
    @test format(-15USD, styles=[:european]) == "−15,00 USD"
    @test format(-1050EUR, styles=[:european]) == "−1.050,00 EUR"
    @test format(5123JPY, styles=[:european]) == "5.123 JPY"
    @test format(12345678.90GBP, styles=[:european]) == "12.345.678,90 GBP"
    @test format(0USD, styles=[:european]) == "0,00 USD"

    # indian style
    @test format(20000000INR, styles=[:indian]) == "INR 2,00,00,000.00"

    # combined us & finance style
    @test format(-1100.55USD, styles=[:us, :finance]) == "USD (1,100.55)"
    @test format(805.03CAD, styles=[:us, :finance]) == "CAD 805.03"
    @test format(1111111.11AUD, styles=[:us, :finance]) == "AUD 1,111,111.11"
    @test format(0USD, styles=[:us, :finance]) == "USD —"

    # combined european & finance styles
    @test format(-1100.55EUR, styles=[:european, :finance]) == "(1.100,55) EUR"
    @test format(-0.11EUR, styles=[:european, :finance]) == "(0,11) EUR"
    @test format(0AUD, styles=[:european, :finance]) == "— AUD"

    # local style
    @test format(1000AUD, styles=[:local]) == "AUD 1000.00"
    @test format(1000EUR, styles=[:local]) == "1000.00 EUR"

    # brief style
    @test format(-15USD, styles=[:brief]) == "−15.00\$"
    @test format(8.05CAD, styles=[:brief]) == "8.05\$"
    @test format(11.11EUR, styles=[:brief]) == "11.11€"

    # combined brief & xxx style
    @test format(-11.11EUR, styles=[:brief, :finance]) == "(11.11€)"
    @test format(-11.11USD, styles=[:brief, :us]) == "−\$11.11"
    @test format(-11.11CAD, styles=[:brief, :european]) == "−11,11\$"
    @test format(-123456.78USD, styles=[:brief, :finance, :us]) ==
        "(\$123,456.78)"
    @test format(805.11GBP, styles=[:brief, :finance, :european]) ==
        "805,11£"

    # combined LaTeX style
    @test format(-11.11AUD, styles=[:latex, :finance]) ==
        "(11.11)\\,\\mathrm{AUD}"
    @test format(11.11AUD, styles=[:latex, :us]) == "\\mathrm{AUD}\\,11.11"
    @test format(8.05AUD, styles=[:latex, :brief, :european]) ==
        "8,05\\mathrm{\\\$}"

    # test declarative nature (duplicates don't matter)
    @testset "Declarative" begin
        for style in (:us, :european, :brief, :finance, :latex)
            for amt in (8000.05AUD, -0.15GBP, 13JPY)
                @test format(amt, styles=[style]) ==
                    format(amt, styles=[style, style])
            end
        end
    end

    # can't combine US & european
    @test_throws Currencies.DeclarativeFormatting.IncompatibleFormatException format(
        USD, styles=[:us, :european])

    # big numbers
    @test format(Monetary{:USD,BigInt,10}(1e4)) == "10000.0000000000 USD"
end


text:
function parse_name!(lexemes::Vector{Token},
                     i::Integer,
                     doc::DocumentationChunk)
    if lexemes[i + 1].kind != WHITESPACE || lexemes[i + 2].kind != TEXT
        print(lexemes[i + 1].contents)
        println(lexemes[i + 2].contents)
        error("Proper format: @@name NAME")
    end

    doc.name = lexemes[i + 2].contents

    if lexemes[i + 3].kind != WHITESPACE || lexemes[i + 4].kind != DIRECTIVE
        print(lexemes[i + 3].contents)
        println(lexemes[i + 4].contents)
        error("Proper format: @@name NAME")
    end

    return i + 4
end

function parse_exported!(lexemes::Vector{Token},
                         i::Integer,
                         doc::DocumentationChunk)
    doc.exported = true

    i += 1

    while i <= length(lexemes) && lexemes[i].kind != DIRECTIVE
        i += 1
    end

    return i
end

function parse_description!(lexemes::Vector{Token},
                            i::Integer,
                            doc::DocumentationChunk)
    if lexemes[i + 1].kind != WHITESPACE ||
       lexemes[i + 2].kind != TEXT
        error("Proper format: @@description DESCRIPTION...")
    end

    io = IOBuffer()

    i = i + 2

    while i <= length(lexemes) && lexemes[i].kind != DIRECTIVE
        if lexemes[i].kind == TEXT
            print(io, lexemes[i].contents)
        else
            print(io, ' ')
        end
        i += 1
    end

    doc.description = bytestring(io)

    return i
end

function parse_arg!(lexemes::Vector{Token},
                      i::Integer,
                      doc::DocumentationChunk)
    if lexemes[i + 1].kind != WHITESPACE ||
       lexemes[i + 2].kind != TEXT ||
       lexemes[i + 3].kind != WHITESPACE
        error("Proper format: @@arg NAME::TYPE DESCRIPTION...")
    end

    varname, vartype = split(lexemes[i + 2].contents, "::")

    io = IOBuffer()

    i = i + 4

    while i <= length(lexemes) && lexemes[i].kind != DIRECTIVE
        if lexemes[i].kind == TEXT
            print(io, lexemes[i].contents)
        else
            print(io, ' ')
        end
        i += 1
    end

    push!(doc.args, Variable(varname, vartype, bytestring(io)))

    return i
end

function parse_field!(lexemes::Vector{Token},
                      i::Integer,
                      doc::DocumentationChunk)
    if lexemes[i + 1].kind != WHITESPACE ||
       lexemes[i + 2].kind != TEXT ||
       lexemes[i + 3].kind != WHITESPACE
        error("Proper format: @@field NAME::TYPE DESCRIPTION...")
    end

    varname, vartype = split(lexemes[i + 2].contents, "::")

    io = IOBuffer()

    i = i + 4

    while i <= length(lexemes) && lexemes[i].kind != DIRECTIVE
        if lexemes[i].kind == TEXT
            print(io, lexemes[i].contents)
        else
            print(io, ' ')
        end
        i += 1
    end

    push!(doc.fields, Variable(varname, vartype, bytestring(io)))

    return i
end

function parse_return!(lexemes::Vector{Token},
                       i::Integer,
                       doc::DocumentationChunk)
    if lexemes[i + 1].kind != WHITESPACE ||
       lexemes[i + 2].kind != TEXT ||
       lexemes[i + 3].kind != WHITESPACE
        error("Proper format: @@return NAME::TYPE DESCRIPTION...")
    end

    varname, vartype = split(lexemes[i + 2].contents, "::")

    io = IOBuffer()

    i = i + 4

    while i <= length(lexemes) && lexemes[i].kind != DIRECTIVE
        if lexemes[i].kind == TEXT
            print(io, lexemes[i].contents)
        else
            print(io, ' ')
        end
        i += 1
    end

    push!(doc.returns, Variable(varname, vartype, bytestring(io)))

    return i
end

function parse_examples!(lexemes::Vector{Token},
                         i::Integer,
                         doc::DocumentationChunk)
    io = IOBuffer()

    i += 1
    while i <= length(lexemes) && lexemes[i].kind != DIRECTIVE
        print(io, lexemes[i].contents)
        i += 1
    end

    doc.examples = bytestring(io)

    return i
end

function parse_directive!(lexemes::Vector{Token},
                          i::Integer,
                          doc::DocumentationChunk)
    parsers = ["@@name" => parse_name!,
               "@@exported" => parse_exported!,
               "@@description" => parse_description!,
               "@@arg" => parse_arg!,
               "@@field" => parse_field!,
               "@@return" => parse_return!,
               "@@examples" => parse_examples!]
    return parsers[lexemes[i].contents](lexemes, i, doc)
end

function parse_plus_semantics(lexemes::Vector{Token})
    chunk = DocumentationChunk()

    # Ignore everything before the first directive
    i = 1
    while lexemes[i].kind != 1 && i <= length(lexemes)
        i += 1
    end

    while i <= length(lexemes)
        if lexemes[i].kind == DIRECTIVE
            i = parse_directive!(lexemes, i, chunk)
        else
            i += 1
        end
    end

    if !isempty(chunk.args) && !isempty(chunk.fields)
        error("Types do not have args; functions do not have fields.")
    end

    return chunk
end


text:
function nestlist( f::Function, x0, n::Integer)
	(n ≤ 0) ? [x0] : begin
		list = Vector{typeof(x0)}(undef, n+1)
		list[1] = x0
		for i in 1:n
			list[i+1] = f(list[i])
		end
		list
	end
end

function bin2dec(bvect::BitVector)
    v = 1
	s = 0
    for i in view(bvect,length(bvect):-1:1)
        s += v*i
        v <<= 1
    end
    s
end


text:
module annulus_examples
using FinEtools
using FinEtoolsHeatDiff
using FinEtoolsHeatDiff.AlgoHeatDiffModule
using LinearAlgebra: cholesky

function annulus_Q4_example()
    println("""
    Annular region,  ingoing and outgoing flux. Minimum/maximum temperature ~(+/-)0.591106.
    Mesh of linear quadrilaterals.
    """)

    t0 = time()

    kappa = 0.2*[1.0 0; 0 1.0]; # conductivity matrix
    magn = 0.06;# heat flux along the boundary
    rin =  1.0;#internal radius
    rex =  2.0;#external radius
    nr = 10; nc = 80;
    Angle = 2*pi;
    thickness =  1.0;
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes  =  Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes  =  mergenodes(fens,  fes,  tolerance);
    edge_fes  =  meshboundary(fes);

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz, 1), 1))

    l1  = selectnode(fens; box=[0.0 0.0 -rex -rex],  inflate = tolerance)
    setebc!(Temp, l1, 1; val=zero(FFlt))
    applyebc!(Temp)

    numberdofs!(Temp)

    material = MatHeatDiff(kappa)
    femm = FEMMHeatDiff(IntegDomain(fes,  GaussRule(2, 2)),  material)

    @time K = conductivity(femm,  geom,  Temp)

    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    el1femm = FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[-magn]);#entering the domain
    @time F1 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    l1 = selectelem(fens, edge_fes, box=[0.9*rex 1.1*rex -0.5*rex 0.5*rex]);
    el1femm =  FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[+magn]);#leaving the domain
    @time F2 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    @time F3 = nzebcloadsconductivity(femm,  geom,  Temp);

    @time K = cholesky(K)
    @time U = K\(F1+F2+F3)
    @time scattersysvec!(Temp, U[:])

    println("Total time elapsed = ", time() - t0, "s")

    File =  "annulus.vtk"
    vtkexportmesh(File,  connasarray(fes),  [geom.values Temp.values], FinEtools.MeshExportModule.Q4; scalars=[("Temperature", Temp.values)])
    @async run(`"paraview.exe" $File`)

    println("Minimum/maximum temperature= $(minimum(Temp.values))/$(maximum(Temp.values)))")

    true

end # annulus_Q4_example

function annulus_Q4_example_algo()
    println("""
    Annular region, ingoing and outgoing flux. Temperature at one node prescribed.
    Minimum/maximum temperature ~(+/-)0.500.
    Mesh of serendipity quadrilaterals.
    This version uses the FinEtools algorithm module.
    Version: 05/29/2017
    """)

    t0 = time()

    kappa = 0.2*[1.0 0; 0 1.0]; # conductivity matrix
    magn = 0.06;# heat flux along the boundary
    rin =  1.0;#internal radius

    rex =  2.0; #external radius
    nr = 3; nc = 40;
    Angle = 2*pi;
    thickness =  1.0;
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes = Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes = mergenodes(fens,  fes,  tolerance);
    edge_fes = meshboundary(fes);

    # At a single point apply an essential boundary condition (pin down the temperature)
    l1  = selectnode(fens; box=[0.0 0.0 -rex -rex],  inflate = tolerance)
    essential1 = FDataDict("node_list"=>l1, "temperature"=>0.0)

    # The flux boundary condition is applied at two pieces of surface
    # Side 1
    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    el1femm = FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[-magn]);#entering the domain
    flux1 = FDataDict("femm"=>el1femm, "normal_flux"=>-magn) # entering the domain
    # Side 2
    l2=selectelem(fens,edge_fes,box=[0.9*rex 1.1*rex -0.5*rex 0.5*rex]);
    el2femm = FEMMBase(IntegDomain(subset(edge_fes, l2),  GaussRule(1, 2)))
    flux2 = FDataDict("femm"=>el2femm, "normal_flux"=>+magn) # leaving the domain

    material = MatHeatDiff(kappa)
    femm = FEMMHeatDiff(IntegDomain(fes,  GaussRule(2, 2)),  material)
    region1 = FDataDict("femm"=>femm)

    # Make model data
    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region1], "essential_bcs"=>[essential1],
    "flux_bcs"=>[flux1, flux2]);

    # Call the solver
    modeldata = AlgoHeatDiffModule.steadystate(modeldata)
    geom=modeldata["geom"]
    Temp=modeldata["temp"]
    println("Minimum/maximum temperature= $(minimum(Temp.values))/$(maximum(Temp.values)))")

    println("Total time elapsed = ",time() - t0,"s")

    # Postprocessing
    vtkexportmesh("annulusmod.vtk", connasarray(fes), [geom.values Temp.values], FinEtools.MeshExportModule.Q4; scalars=[("Temperature", Temp.values)])

end # annulus_Q4_example_algo

function annulus_Q8_example()
    println("""
    Annular region,  ingoing and outgoing flux. Minimum/maximum temperature ~(+/-)0.501.
    Mesh of quadratic serendipity quadrilaterals.
    Version: 05/29/2017
    """)

    t0 = time()

    kappa = 0.2*[1.0 0; 0 1.0]; # conductivity matrix
    magn = 0.06;# heat flux along the boundary
    rin =  1.0;#internal radius

    rex =  2.0; #external radius
    nr = 5; nc = 40;
    Angle = 2*pi;
    thickness =  1.0;
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes  =  Q8annulus(rin, rex, nr, nc, Angle)
    fens, fes  =  mergenodes(fens,  fes,  tolerance);
    edge_fes  =  meshboundary(fes);

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz, 1), 1))

    l1  = selectnode(fens; box=[0.0 0.0 -rex -rex],  inflate = tolerance)
    setebc!(Temp, l1, 1; val=zero(FFlt))
    applyebc!(Temp)

    numberdofs!(Temp)

    material = MatHeatDiff(kappa)
    femm = FEMMHeatDiff(IntegDomain(fes,  GaussRule(2, 3)),  material)

    @time K = conductivity(femm,  geom,  Temp)

    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    el1femm = FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[-magn]);#entering the domain
    @time F1 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    l1 = selectelem(fens, edge_fes, box=[0.9*rex 1.1*rex -0.5*rex 0.5*rex]);
    el1femm =  FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[+magn]);#leaving the domain
    @time F2 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    @time F3 = nzebcloadsconductivity(femm,  geom,  Temp);

    @time K = cholesky(K)
    @time U = K\(F1+F2+F3)
    @time scattersysvec!(Temp, U[:])

    println("Total time elapsed = ", time() - t0, "s")

    File =  "annulusq8.vtk"
    vtkexportmesh(File, connasarray(fes),  [geom.values Temp.values], FinEtools.MeshExportModule.Q8; scalars=[("Temperature", Temp.values)])

    println("Minimum/maximum temperature= $(minimum(Temp.values))/$(maximum(Temp.values)))")

    true

end # annulus_Q8_example

function ebc_annulus_Q4_algo()
    println("""
    Annular region, ingoing and outgoing flux. Temperature at one node prescribed.
    Minimum/maximum temperature ~(+/-)0.500.
    Mesh of serendipity quadrilaterals.
    This version uses the FinEtools algorithm module.
    Version: 05/29/2017
    """)

    t0 = time()

    kappa = 0.2*[1.0 0; 0 1.0]; # conductivity matrix
    te1 = -0.5
    te2 = 0.6
    hconv1, hconv2 = 1000.0, 1000.0
    rin =  1.0;#internal radius
    rex =  2.0; #external radius
    nr = 7; nc = 90;
    Angle = 2*pi;
    thickness =  1.0;
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes = Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes = mergenodes(fens,  fes,  tolerance);
    edge_fes = meshboundary(fes);

    # The convection boundary condition is applied at two pieces of surface
    # Side 1
    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    el1femm = FEMMHeatDiffSurf(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)), hconv1)
    cbc1 = FDataDict("femm"=>el1femm, "ambient_temperature"=>te1)
    # Side 2
    l2=selectelem(fens,edge_fes,box=[0.9*rex 1.1*rex -0.5*rex 0.5*rex]);
    el2femm = FEMMHeatDiffSurf(IntegDomain(subset(edge_fes, l2),  GaussRule(1, 2)), hconv2)
    cbc2 = FDataDict("femm"=>el2femm, "ambient_temperature"=>te2)

    material = MatHeatDiff(kappa)
    femm = FEMMHeatDiff(IntegDomain(fes,  GaussRule(2, 2)),  material)
    region1 = FDataDict("femm"=>femm)

    # Make model data
    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region1],
    "convection_bcs"=>[cbc1, cbc2]);

    # Call the solver
    modeldata = AlgoHeatDiffModule.steadystate(modeldata)
    geom=modeldata["geom"]
    Temp=modeldata["temp"]
    println("Minimum/maximum temperature= $(minimum(Temp.values))/$(maximum(Temp.values)))")

    println("Total time elapsed = ",time() - t0,"s")

    # Postprocessing
    File = "annulusmod.vtk"
    vtkexportmesh(File, connasarray(fes), [geom.values Temp.values],  FinEtools.MeshExportModule.Q4; scalars=[("Temperature", Temp.values)])
    @async run(`"paraview.exe" $File`)

end # ebc_annulus_Q4_algo

function allrun()
    println("#####################################################")
    println("# annulus_Q4_example ")
    annulus_Q4_example()
    println("#####################################################")
    println("# annulus_Q4_example_algo ")
    annulus_Q4_example_algo()
    println("#####################################################")
    println("# annulus_Q8_example ")
    annulus_Q8_example()
    println("#####################################################")
    println("# ebc_annulus_Q4_algo ")
    ebc_annulus_Q4_algo()
end # function allrun

end # module annulus_examples

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
#!/usr/bin/env julia

using JLD2
using FileIO
using SymmetricTensors
using PyCall
@pyimport matplotlib as mpl
using PyPlot
mpl.rc("text", usetex=true)
mpl.rc("font", family="serif", size = 8)

"""
  pltdiag()

Plots a chart of superdiagonal elements of cumulants of and its theoretical values
"""
function pltdiag()
  tdiag = try load("datafortests.jld2")["theoretical diag"]
  catch
    println("please run test/gandata.jl and test/testondata.jl")
  return ()
  end
  cum = try load("cumulants.jld2")["cumulants"]
  catch
    println("please run test/testondata.jl")
    return ()
  end
  n = cum[1].dats
  fig, ax = subplots(figsize = (3., 2.3))
  col = ["cyan", "brown", "green", "red", "blue", "black"]
  for order in (2,4,5,6)
    c = cum[order]
    ax[:plot](diag(c), "o", color = col[order], label = "$order cumulant", markersize=3)
    ax[:plot]([fill(tdiag[order], n)...], "--", color = col[order], label = "theoretical")
  end
  PyPlot.ylabel("superdiagonal elements", labelpad = -1)
  PyPlot.xlabel("superdiagonal index", labelpad = -3)
  ax[:legend](fontsize = 4.5, loc = 5)
  fig[:savefig]("diagcumels.pdf")
end

function main()
  pltdiag()
end

main()


text:
using Test
using DojoEnvironments

using Dojo

@testset "DojoEnvironments" begin
    @testset "mechanisms" begin
        include("mechanisms.jl")
    end
    @testset "environments" begin
        include("environments.jl")
    end
end


text:
"""
    EntropyBasedCausalityTest

The supertype of all abstract and composite types representing a causality
test based on some entropy based measure.

Concrete subtypes are those based on transfer entropy

- [`VisitationFrequencyTest`](@ref)
- [`TransferOperatorGridTest`](@ref)
- [`ExactSimplexIntersectionTest`](@ref)
- [`ApproximateSimplexIntersectionTest`](@ref)

and those that in some manner utilise transfer entropy or other
information theoretic approaches :

- [`PredictiveAsymmetryTest`](@ref)
"""
abstract type EntropyBasedCausalityTest <: CausalityTest end

export EntropyBasedCausalityTest


text:
function Base.sum(r::StaticStepRangeLen)
    l = length(r)
    # Compute the contribution of step over all indices.
    # Indexes on opposite side of r.offset contribute with opposite sign,
    #    r.step * (sum(1:np) - sum(1:nn))
    np, nn = l - offset(r), offset(r) - SOne  # positive, negative
    # To prevent overflow in sum(1:n), multiply its factors by the step
    sp, sn = Base.sumpair(np), Base.sumpair(nn)
    tp = _tp_prod(step(r), sp[1], sp[2])
    tn = _tp_prod(step(r), sn[1], sn[2])
    s_hi, s_lo = add12(gethi(tp), -gethi(tn))
    s_lo += getlo(tp) - getlo(tn)
    # Add in contributions of ref
    ref = reference(r) * l
    sm_hi, sm_lo = add12(s_hi, gethi(ref))
    add12(sm_hi, sm_lo + getlo(ref))[1]
end

function Base._tp_prod(t::TPVal, x, y...)
    Base.@_inline_meta
    Base._tp_prod(t * x, y...)
end
Base._tp_prod(t::TPVal) = t

function +(r1::StaticStepRangeLen{T,R1}, r2::StaticStepRangeLen{T,R2}
          ) where {T,R1<:Union{<:TPVal,TwicePrecision},R2<:Union{<:TPVal,TwicePrecision}}
    len = length(r1)
    (len == length(r2) || throw(DimensionMismatch("argument dimensions must match")))
    if offset(r1) == offset(r2)
        imid = offset(r1)
        ref = reference(r1) + reference(r2)
    else
        imid = round(Int, (offset(r1)+offset(r2))/ofeltype(2))
        ref1mid = _getindex_hiprec(r1, imid)
        ref2mid = _getindex_hiprec(r2, imid)
        ref = ref1mid + ref2mid
    end
    step = twiceprecision(step(r1) + step(r2), nbitslen(T, len, imid))
    StaticStepRangeLen{T,typeof(ref),typeof(step)}(ref, step, len, imid)
end

Base.sumpair(n::SInteger) = iseven(n) ? (n+SOne, n>>SOne) : (n, (n+SOne)>>SOne)


text:
@testset "1631.path-with-minimum-effort.jl" begin
    @test minimum_effort_path([[1,2,2],[3,8,2],[5,3,5]]) == 2
    @test minimum_effort_path([[1,2,3],[3,8,4],[5,3,5]]) == 1
    @test minimum_effort_path([[1,2,1,1,1],[1,2,1,2,1],[1,2,1,2,1],[1,2,1,2,1],[1,1,1,2,1]]) == 0
    @test minimum_effort_path([[1,10,6,7,9,10,4,9]]) == 9
end