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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ d7bf8b8d-de6b-4783-895c-e2f5a1850067
begin
	import Pkg
	Pkg.activate(".")
	# Pkg.add("Plots")
	# Pkg.add("CSV")
	# Pkg.add("DataFrames")
	# Pkg.add("Dates")
	# Pkg.add("ImageFiltering")
	# Pkg.add("OffsetArrays")
	# Pkg.add("Statistics")
	# Pkg.add("PlutoUI")
end

# ╔═╡ ef23d264-448e-11eb-10f8-ad25ae5f244d
begin
	using Plots
	using CSV
	using DataFrames
	using Dates
	using ImageFiltering
	using OffsetArrays
	using Statistics
	using PlutoUI
end

# ╔═╡ e2635196-448c-11eb-2877-c5735be35d80
md"""
## California COVID correlations

This notebook will take COVID case, hospitalization, and death data from California counties and attempt to find a lagged CFR in the style of Trevor Bedford's [analysis](https://twitter.com/trvrb/status/1326404889438801921).

The main idea is that a reasonably stable fraction of cases on a given day will die from COVID but that it takes some time for the disease to progress, so the death rate lags the case rate. Of course, we don't assume that we are identifying all the cases, so the lagged case-fatality-rate (CFR) here is not a real CFR. It is likely a big overestimate, but if it is relatively stable, it will allow us to forecast deaths 15+ days ahead of today.
"""

# ╔═╡ cdec3534-43f1-4726-af95-d82e32f3a0cb
begin
	state_case_data = CSV.read("covid19cases_test.csv", DataFrame)
	county_list = unique(state_case_data[!, :area])
	state_hosp_data = CSV.read("statewide-covid-19-hospital-county-data.csv", DataFrame)
	county_list_2 = unique(state_hosp_data[!, :county])
	county_list = intersect(county_list, county_list_2)
end;

# ╔═╡ 9887d642-8387-477c-991c-4756bc137252
md"""
Using the dropdown, select the county to investigate. This analysis works best with big populations, so smaller counties may not work as well.
"""

# ╔═╡ 3ac6e767-2df5-4969-8582-11312228f534
md"""
## Data from $(@bind county_name Select(county_list, default="Santa Clara")) county

"""

# ╔═╡ 1d7590fd-8295-4964-b368-4629b1b3137c
begin
	offset = 3 #recent data is somewhat unreliable
	county_case_data  = filter(row -> row.area == county_name, state_case_data)
	county_case_dates = Date.(county_case_data[!, "date"][1:end-offset])
	county_cases = county_case_data[!, "cases"][1:end-offset]
	county_deaths = county_case_data[!, "deaths"][1:end-offset]
end;

# ╔═╡ 7a75c70f-e7d5-4505-8b6e-128e3a9c2e0b
begin
	county_hosp_data  = filter(row -> row.county == county_name, state_hosp_data)
	county_hosp_dates = Date.(county_hosp_data[!, "todays_date"])
	county_tot_patients= (county_hosp_data[!, "hospitalized_covid_confirmed_patients"] 							)#.+county_hosp_data[!, "hospitalized_suspected_covid_patients"])

	county_icu_cases = (county_hosp_data[!, "icu_covid_confirmed_patients"] 							.+county_hosp_data[!, "icu_suspected_covid_patients"])
	county_non_icu = county_tot_patients-county_icu_cases
end;

# ╔═╡ 8f0fce44-46ed-11eb-395c-b3072d1362e3
md"""
As one can see below, the case reporting data has artificial feautures due to a relative lack of reporting on Sundays. To correct for this, we'll take the 7-day average looking backwards.
"""

# ╔═╡ 9685d04e-4493-11eb-0db5-c35b82182d7f
begin
	StartDate = Date(2020,6,1);
	case_start_ind = findfirst(x -> x==StartDate, county_case_dates);
	hosp_start_ind = findfirst(x -> x==StartDate, county_hosp_dates);
end;

# ╔═╡ f7923db6-448f-11eb-0e94-e77b1a911cff
begin
	plotly()
	plot(county_case_dates[case_start_ind:end],
		county_cases[case_start_ind:end], label = "New Cases")
end

# ╔═╡ 2fac8934-4494-11eb-0219-7f1d44649362
kernel = OffsetArray(fill(1/7, 7), -6:0);

# ╔═╡ 110729c7-6d6e-4efc-939f-d945a6d5ee12
county_moving_avg = imfilter(county_cases[case_start_ind:end], kernel);

# ╔═╡ 257048ea-4495-11eb-0b21-f353c7997010
begin
	plotly()
	plot(county_case_dates[case_start_ind:end], county_moving_avg, label = "Avg. Cases", legend=:topleft)
	plot!(county_case_dates[case_start_ind:end],
		county_cases[case_start_ind:end], label = "New Cases")
end

# ╔═╡ 4efe6d30-4498-11eb-2e05-6924e9b3aa6a
begin
	county_avg_deaths = imfilter(county_deaths[case_start_ind:end],kernel)
	compMovingAvg = imfilter(county_cases[case_start_ind:end], kernel);
end;

# ╔═╡ 24d95b88-0521-49b6-9258-a96b10a8e541
md"""
## Comparing lagged cases with deaths
"""

# ╔═╡ f6f679d8-449a-11eb-1967-3da2b02ef0dc
@bind lag Slider(5:30, default=17, show_value=true)

# ╔═╡ 681f0498-449b-11eb-2dee-7dd6d54dfc4e
@bind lag_cfr Slider(0.001:0.001:0.03, default=0.018, show_value=true)

# ╔═╡ eed3fb32-46ed-11eb-13f6-0d79c4527fb0
md"""
Given the population of $(county_name) county, the average deaths are still fairly jagged because we still have small numbers. This death count in the county itself doesn't track the lagged case rate as well as people have seen using larger populations. Here I'm guessing we're seeing the inhomogeneity of COVID, e.g., cases counts could go up in low-risk populations and thus not see deaths rise until quite a bit later or conversely a small outbreak at a long-term care facility could lead to a lot of deaths.

Still the general trend is there, when cases rise sharply, deaths rise about $lag days later. The current fit predicts that for 100 new, identified cases, $(round((lag_cfr*100)*10)/10) people will die from COVID.
"""

# ╔═╡ 5b04aebe-4498-11eb-1309-93cbca18c737
begin
	plot(county_case_dates[case_start_ind:end]+Day(lag), lag_cfr*compMovingAvg[1:end],
		label = "Estimated deaths", legend=:topleft)
	plot!(county_case_dates[case_start_ind:end],
		county_avg_deaths, label = "County deaths")
end

# ╔═╡ a14162e8-46ee-11eb-3edb-8f55a7565c4b
md"""
Having made our rough fit, we can project deaths going foreward by taking recent case counts, multiplying by the lagged CFR and shifting forward in time. Below we see that by $(county_case_dates[end-5]+Day(lag)), there should be roughly $(round(lag_cfr*compMovingAvg[end-5])) deaths per day.

This prediction should be an overestimate given the presence of vaccines, but if the cases are concentrated in the unvaccinated, it might not be.
"""

# ╔═╡ 35eb6a82-449a-11eb-1f5e-ebe3e5fa4fb6
begin
	plotly()
	plot(county_case_dates[end-lag:end-5]+Day(lag),
		lag_cfr*compMovingAvg[end-lag:end-5],
		label = "Estimated deaths",
		legend=:topleft)
end

# ╔═╡ 291c556e-44b6-11eb-0666-df947594a060
md"""
## Comparing lagged cases to hospitalization
Now we're going to add in the number of hospitalized people. By similar arguments, a fraction of cases on a given day will eventually need hospital beds or potentially an ICU bed. We can also fit the averaged new case data to these figures as well.

Note, the adjustment factor here is not a simple hospitalization rate since people stay for multiple days in the ICU or hospital. If we knew the average time a COVID patient spends in the ICU or hospital, we could figure out the case/hospital bed or case/ICU ratio, but this would still have the issues from lack of testing and thus could not be used to predict how frequently a COVID infection requires hospital services.

Nevertheless, we can still use these correlation to predict future hospital need.
"""

# ╔═╡ ba472a38-46ef-11eb-21e6-97331ae5e715
md"""
Below, we see much better agreement between the scaled, lagged case numbers and the number of ICU patients.

However, somewhere in mid-December, just as ICUs began to approach capacity, we can see a clear change in the rate of people admitted. This may have arisen for many reasons, two of the most obvious would be only admitting the most severe cases to the ICU to reserve beds or an increase in cases from increased testing as people prepared for the holidays. The degree to which the deaths tracked cases during the same period could provide some evidence for the former vs. the latter.
"""

# ╔═╡ 3fd3dce4-44b8-11eb-0f23-93765a452257
@bind icu_lag Slider(1:20, default=15, show_value=true)

# ╔═╡ f8fffc76-44b7-11eb-3de5-9d2a8a97ca11
@bind icu_rate Slider(0.0:0.01:.3, default = 0.19, show_value=true)

# ╔═╡ 2a54ec88-44b7-11eb-196d-f7844a62e8f6
begin
	plotly()
	plot(county_hosp_dates[hosp_start_ind:end], county_icu_cases[hosp_start_ind:end], label = "ICU patients", legend=:topleft)
	plot!(county_case_dates[case_start_ind:end]+Day(icu_lag),
		icu_rate*compMovingAvg,
		label = "Scaled, lagged cases", legend=:topleft)
	#plot!(county_dates, county_icu_cases)
end

# ╔═╡ b5fb79f2-44bf-11eb-025b-ed9492e39ade
@bind hosp_lag Slider(1:20, default=13, show_value=true)

# ╔═╡ 30323758-44c0-11eb-0df3-1ddd46872b74
md"""
If we assume the average ICU patient spends 5 days there, we predict that ~$(round(icu_rate*100/5)) % of identified COVID cases need ICU beds.

We can also predict that on $(county_case_dates[end]+Day(hosp_lag)) there will be $(round(icu_rate*compMovingAvg[end])) ICU patients hospitalized with COVID19
"""

# ╔═╡ 03f58904-46f0-11eb-08bb-035012ae329d
md"""
### Non-ICU patients
We repeat the analysis with non-ICU COVID patients and see our best correlation yet. The lag is only $hosp_lag days, so our predictive powers are limited, however they show a clear indication that rising cases leads to rising COVID+ patients who are hospitalized.
"""

# ╔═╡ c0b2387c-44bf-11eb-31d1-3f403a6067f3
@bind hosp_rate Slider(0.0:0.01:.6, default = 0.47, show_value=true)

# ╔═╡ f37730a2-44b8-11eb-18d2-670aa5a0209e
begin
	plot(county_hosp_dates[hosp_start_ind:end],
		county_non_icu[hosp_start_ind:end], label= "Non-ICU Patients")
	plot!(county_case_dates[case_start_ind:end]+Day(hosp_lag),
		hosp_rate*compMovingAvg[1:end],
		label = "Scaled, lagged cases", legend=:topleft)

end

# ╔═╡ f3525c1c-44bf-11eb-0619-4f99561f4a6b
md"""
If we assume the average hospital stay is 5 days, each identified COVID case as a $(round(hosp_rate/5*100)) % chance of needing a hospital bed.
"""

# ╔═╡ ede8c6f0-51d8-11eb-2681-cb27cee40af5
md"""
Current prediction is that on $(county_case_dates[end]+Day(hosp_lag)) there will be $(round(hosp_rate*compMovingAvg[end])) non-ICU patients hospitalized with COVID19
"""

# ╔═╡ b79d0f23-c592-4e03-b9e9-dc61b3cedb6e
md"""
#### Import statements and helper code
"""

# ╔═╡ 2f74be36-44ae-11eb-278d-9db0cdf4a98b
begin
	state_cases_deaths_url = "https://data.chhs.ca.gov/dataset/f333528b-4d38-4814-bebb-12db1f10f535/resource/046cdd2b-31e5-4d34-9ed3-b48cdbc4be7a/download/covid19cases_test.csv"
	download(state_cases_deaths_url,  "covid19cases_test.csv")
	state_hosp_url = "https://data.chhs.ca.gov/dataset/2df3e19e-9ee4-42a6-a087-9761f82033f6/resource/47af979d-8685-4981-bced-96a6b79d3ed5/download/covid19hospitalbycounty.csv"
	download(state_hosp_url,  "covid19hospitalbycounty.csv")
end;

# ╔═╡ Cell order:
# ╟─e2635196-448c-11eb-2877-c5735be35d80
# ╠═d7bf8b8d-de6b-4783-895c-e2f5a1850067
# ╠═ef23d264-448e-11eb-10f8-ad25ae5f244d
# ╠═cdec3534-43f1-4726-af95-d82e32f3a0cb
# ╟─9887d642-8387-477c-991c-4756bc137252
# ╟─3ac6e767-2df5-4969-8582-11312228f534
# ╠═1d7590fd-8295-4964-b368-4629b1b3137c
# ╠═7a75c70f-e7d5-4505-8b6e-128e3a9c2e0b
# ╟─8f0fce44-46ed-11eb-395c-b3072d1362e3
# ╟─f7923db6-448f-11eb-0e94-e77b1a911cff
# ╠═9685d04e-4493-11eb-0db5-c35b82182d7f
# ╠═2fac8934-4494-11eb-0219-7f1d44649362
# ╠═110729c7-6d6e-4efc-939f-d945a6d5ee12
# ╠═257048ea-4495-11eb-0b21-f353c7997010
# ╟─4efe6d30-4498-11eb-2e05-6924e9b3aa6a
# ╟─24d95b88-0521-49b6-9258-a96b10a8e541
# ╟─eed3fb32-46ed-11eb-13f6-0d79c4527fb0
# ╠═f6f679d8-449a-11eb-1967-3da2b02ef0dc
# ╠═681f0498-449b-11eb-2dee-7dd6d54dfc4e
# ╟─5b04aebe-4498-11eb-1309-93cbca18c737
# ╟─a14162e8-46ee-11eb-3edb-8f55a7565c4b
# ╟─35eb6a82-449a-11eb-1f5e-ebe3e5fa4fb6
# ╟─291c556e-44b6-11eb-0666-df947594a060
# ╟─ba472a38-46ef-11eb-21e6-97331ae5e715
# ╠═3fd3dce4-44b8-11eb-0f23-93765a452257
# ╠═f8fffc76-44b7-11eb-3de5-9d2a8a97ca11
# ╠═2a54ec88-44b7-11eb-196d-f7844a62e8f6
# ╟─30323758-44c0-11eb-0df3-1ddd46872b74
# ╟─03f58904-46f0-11eb-08bb-035012ae329d
# ╠═b5fb79f2-44bf-11eb-025b-ed9492e39ade
# ╠═c0b2387c-44bf-11eb-31d1-3f403a6067f3
# ╟─f37730a2-44b8-11eb-18d2-670aa5a0209e
# ╟─f3525c1c-44bf-11eb-0619-4f99561f4a6b
# ╟─ede8c6f0-51d8-11eb-2681-cb27cee40af5
# ╟─b79d0f23-c592-4e03-b9e9-dc61b3cedb6e
# ╠═2f74be36-44ae-11eb-278d-9db0cdf4a98b

max_stars_repo_name: sragole/COVIDCorrelations

max_issues_repo_name: sragole/COVIDCorrelations

max_forks_repo_name: sragole/COVIDCorrelations


text:
@testset "Test Evaluator" begin
    @testset "General" begin
        # General tests for all backends
        test_backend(m.evaluate, "evaluator")
    end

    # These cases should not occur unless the internal data are malformed.
    @testset "Test Malformed Expressions" begin
        test_object(
            m.evaluate_prefix_expression("+", m.IntegerObj(1)),
            "unknown operator: +INTEGER",
        )

        test_object(
            m.evaluate_infix_expression("&", m.IntegerObj(1), m.IntegerObj(1)),
            "unknown operator: INTEGER & INTEGER",
        )
    end

    @testset "Test Defining Macro" begin
        code = """
          let number = 1;
          let function = fn(x, y) { x + y };
          let mymacro = macro(x, y) { x + y };
        """

        env = m.Environment()
        program = m.define_macros!(env, m.parse(code))

        @test length(program.statements) == 2
        @test isnothing(m.get(env, "number"))
        @test isnothing(m.get(env, "function"))

        mc = m.get(env, "mymacro")
        @test !isnothing(mc)
        @test length(mc.parameters) == 2
        @test string(mc.parameters[1]) == "x"
        @test string(mc.parameters[2]) == "y"
        @test string(mc.body) == "(x + y)"
    end

    @testset "Test Expanding Macros" begin
        for (code, expected) in [
            (
                "let infixExpression = macro() { quote(1 + 2) }; infixExpression();",
                "(1 + 2)",
            ),
            (
                "let reverse = macro(a, b) { quote(unquote(b) - unquote(a)); }; reverse(2 + 2, 10 - 5);",
                "((10 - 5) - (2 + 2))",
            ),
            (
                """
                let unless = macro(condition, consequence, alternative) {
                    quote(if (!(unquote(condition))) {
                        unquote(consequence);
                    } else {
                        unquote(alternative);
                    });
                };

                unless(10 > 5, puts("not greater"), puts("greater"));
                """,
                "if ((!(10 > 5))) { puts(\"not greater\") } else { puts(\"greater\") }",
            ),
        ]
            env = m.Environment()
            program = m.define_macros!(env, m.parse(code))
            expanded = m.expand_macros(program, env)
            @test string(expanded) == expected
        end
    end

    @testset "String macro" begin
        for (code, expected, expected_output) in [
            ("let b = 4; b;", 4, ""),
            ("puts([1, 2, 3]);", nothing, "[1, 2, 3]\n"),
            ("let m = macro(x, y) { quote(unquote(y) - unquote(x)) }; m(5, 10);", 5, ""),
        ]
            c = IOCapture.capture() do
                eval(quote
                    m.@monkey_eval_str($code)
                end)
            end

            test_object(c.value, expected)
            @test c.output == expected_output
        end
    end
end

ext: jl


text:
# 2D FWI gradient test with 1 source
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

parsed_args = parse_commandline()

nlayer = parsed_args["nlayer"]
tti = parsed_args["tti"]
fs =  parsed_args["fs"]

### Model
model, model0, dm = setup_model(parsed_args["tti"], 4)
q, srcGeometry, recGeometry, info = setup_geom(model)
dt = srcGeometry.dt[1]

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

@testset "FWI gradient test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
	# Gradient test
	h = 5f-2
	maxiter = 6
	err1 = zeros(maxiter)
	err2 = zeros(maxiter)
	h_all = zeros(maxiter)
	modelH = deepcopy(model0)

	# Observed data
	opt = Options(sum_padding=true, free_surface=parsed_args["fs"])
	F = judiModeling(info, model, srcGeometry, recGeometry; options=opt)
	d = F*q

	# FWI gradient and function value for m0
	Jm0, grad = fwi_objective(model0, q, d; options=opt)

	dJ = dot(grad, dm)

	for j=1:maxiter
		# FWI gradient and function falue for m0 + h*dm
		modelH.m = model0.m + h*dm
		Jm, gradm = fwi_objective(modelH, q, d;options=opt)

		# Check convergence
		err1[j] = abs(Jm - Jm0)
		err2[j] = abs(Jm - (Jm0 + h*dJ))
		j == 1 ? prev = 1 : prev = j - 1
		@printf("h = %2.2e, e1 = %2.2e, rate = %2.2e", h, err1[j], err1[prev]/err1[j])
		@printf(", e2  = %2.2e, rate = %2.2e \n", err2[j], err2[prev]/err2[j])
		h_all[j] = h
		h = h * .8f0
	end

	# CHeck convergence rates
	rate_1 = sum(err1[1:end-1]./err1[2:end])/(maxiter - 1)
	rate_2 = sum(err2[1:end-1]./err2[2:end])/(maxiter - 1)

	@test isapprox(rate_1, 1.25f0; rtol=5f-2)
	@test isapprox(rate_2, 1.5625f0; rtol=5f-2)

end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using FractionalCalculus
using Test

@testset "Test Caputo Direct Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 0, 1, 1e-8, Caputo_Direct())[1], 2/sqrt(pi); atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 0, 3.2, 1e-8, Caputo_Direct())[1], 4.300306216488329e2; atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.25, 0, 3.2, 1e-8, Caputo_Direct())[1], 382.1228113951680; atol = 1e-5)
    @test isapprox(fracdiff(x->(x-3)^3+(x+2)^2+5, 0.5, 0, 0.83, 1e-8, Caputo_Direct())[1], 23.899930581118699; atol = 1e-5)
    @test isapprox(fracdiff(x->(x-3)^3+(x+2)^2+5, 0.25, 0, 0.83, 1e-8, Caputo_Direct())[1], 22.96354626572943; atol = 1e-5)
    @test isapprox(fracdiff(x->cos(x), 0.5, 0, 0.34, 1e-8, Caputo_Direct())[1], -0.147174774607683; atol = 1e-5)
    @test isapprox(fracdiff(x->cos(x), 0.25, 0, 0.34, 1e-8, Caputo_Direct())[1], -0.093074344902431; atol = 1e-5)
    @test isapprox(fracdiff(x->exp(-x^2), 0.5, 0, 2.3, 1e-8, Caputo_Direct())[1], -0.505891017154289; atol = 1e-5)
    @test isapprox(fracdiff(x->exp(-x^2), 0.25, 0, 2.3, 1e-8, Caputo_Direct())[1], -0.762927553656252; atol = 1e-5)
end

@testset "Test Caputo Piecewise Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 1e-4, Caputo_Piecewise()), 2/sqrt(pi); atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 3.2, 1e-4, Caputo_Piecewise()), 4.300306216488329e2; atol = 1e-5)
    @test isapprox(fracdiff(cos, 0.5, 0.34, 1e-4, Caputo_Piecewise()), -0.147174774607683; atol = 1e-5)
    @test isapprox(fracdiff(x->exp(-x^2), 0.25, 2.3, 1e-3, Caputo_Piecewise()), -0.762927553656252; atol = 1e-5)
end

@testset "Test Caputo Diethelm Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 1e-4, Caputo_Diethelm()), 2/sqrt(pi); atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 3.2, 1e-5, Caputo_Diethelm()), 4.300306216488329e2; atol = 1e-5)
end

@testset "Test Grunwald-Letnikov Fractional Derivative" begin
    @test isapprox(fracdiff(x->x, 0.5, 0, 0.5, GL_Direct())[1], 0.7978845583186518; atol = 1e-5)
    @test isapprox(fracdiff(x->x^5, 0.5, 0, 3.2, GL_Direct())[1], 4.300306216488329e2; atol = 1e-5)
    @test isapprox(fracdiff(x->x, 0.5, collect(0:0.01:1), 2, GL_High_Precision())[end], 2/sqrt(pi); atol = 1e-4)
end

@testset "Test RLDiff Matrix Fractional Deriavtive" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.001, RLDiff_Matrix())[end], 2/sqrt(pi); atol = 1e-3)
end

@testset "Test RLDiff Approx" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.001, RLDiff_Approx()), 2/sqrt(pi); atol = 1e-4)

    # Test Array input end point
    result = fracdiff(x->x, 0.5, [1, 2, 3], 0.001, RLDiff_Approx())
    @test isapprox(result[1], 2/sqrt(pi); atol=1e-4)
    @test isapprox(result[2], 1.5957691216057306; atol=1e-4)
    @test isapprox(result[3], 1.954410047611678; atol=1e-4)
end

@testset "Test RL Linear Spline Interpolation" begin
    @test isapprox(fracdiff(x->x, 0.5, 1, 0.0001, RL_Linear_Spline_Interp()), 2/sqrt(pi); atol = 1e-4)

    # Test Array input end point
    result = fracdiff(x->x, 0.5, [1, 2, 3], 0.001, RL_Linear_Spline_Interp())
    @test isapprox(result[1], 2/sqrt(pi); atol=1e-4)
    @test isapprox(result[2], 1.5957691216057306; atol=1e-4)
    @test isapprox(result[3], 1.954410047611678; atol=1e-4)
end

@testset "Test Hadamard Fractional Derivative" begin
    @test isapprox(fracdiff(log, 0.5, 1, 2, 1/4000000, Hadamard_LRect()), 0.9391460; atol=1e-4)
    @test isapprox(fracdiff(log, 0.3, 1, 2, 1/4000000, Hadamard_LRect()), 0.8514935; atol=1e-4)
    @test isapprox(fracdiff(log, 0.5, 1, 2, 1/4000000, Hadamard_RRect()), 0.9391460; atol=1e-4)
    @test isapprox(fracdiff(log, 0.3, 1, 2, 1/4000000, Hadamard_RRect()), 0.8514936; atol=1e-4)
    @test isapprox(fracdiff(log, 0.5, 1, 2, 1/4000000, Hadamard_Trap()), 0.9391460; atol=1e-4)
    @test isapprox(fracdiff(log, 0.3, 1, 2, 1/4000000, Hadamard_Trap()), 0.8514935; atol=1e-4)
end

@testset "Test macros" begin
    @test isapprox(@fracdiff(x->x, 0.5, 1), 1.1283791670955126; atol=1e-4)
    @test isapprox(@fracdiff(x->x^5, 0.5, 3.2), 4.300306216488329e2; atol=1e-2)
end


text:
using JuMP
using CPLEX
include("parser.jl")

function dualisation(MyFileName::String)
  n, s, t, S, d1, d2, p, ph, d, D = read_instance(MyFileName)

  # Create the model
  m = Model(CPLEX.Optimizer)

  ## Variables
  @variable(m, x[1:n,1:n], Bin) # certains arcs n'existent pas, ils ne sont pas utilises dans les contraintes
  @variable(m, y[1:n], Bin)
  @variable(m, t1 >= 0)
  @variable(m, t2 >= 0)
  @variable(m, z[1:n,1:n] >= 0, Int) # meme remarque
  @variable(m, zprim[1:n] >= 0, Int)

  ## Constraints
  @constraint(m,sum(x[s,j] for j in 1:n if d[s,j] != 0) == 1)
  @constraint(m,sum(x[i,t] for i in 1:n if d[i,t] != 0) == 1)
  @constraint(m,[v in 1:n; v != s && v != t], sum(x[i,v] for i in 1:n if d[i,v] != 0) == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m,[v in 1:n; v != t], y[v] == sum(x[v,j] for j in 1:n if d[v,j] != 0))
  @constraint(m,[i in 1:n; i != t], y[t] == sum(x[i,t] for i in 1:n if d[i,t] != 0))
  @constraint(m, t2*d2 + sum(p[i]*y[i] + 2*zprim[i] for i in 1:n) <= S)
  @constraint(m, [i in 1:n, j in 1:n ; d[i,j] != 0], t1 + z[i,j] >= d[i,j]*x[i,j])
  @constraint(m, [i in 1:n], t2 + zprim[i] >= ph[i]*y[i])

  ## Objective
  @objective(m, Min, d1*t1 + sum(d[i,j]*x[i,j] + D[i,j]*z[i,j] for i in 1:n, j in 1:n if d[i,j] != 0))

  #resolution
  optimize!(m)

  object = JuMP.objective_value.(m)
  x_star = JuMP.getvalue.( m[:x] )
  y_star = JuMP.getvalue.( m[:y] )
  t1_star = JuMP.getvalue.( m[:t1] )
  t2_star = JuMP.getvalue.( m[:t2] )
  z_star = JuMP.getvalue.( m[:z] )
  zprim_star = JuMP.getvalue.( m[:zprim] )
  println("objectif ", object)
  println("x_star ", x_star)
  println("y_star ", y_star)
  println("z_star ", z_star)
  println("t1_star ", t1_star)
  println("t2_star ", t2_star)
  println("zprim_star ", zprim_star)

end

max_stars_repo_name: boumpteryx/ECMA

max_issues_repo_name: boumpteryx/ECMA

max_forks_repo_name: boumpteryx/ECMA


text:
abstract type AbstractKKTSystem{T, MT<:AbstractMatrix{T}} end

"Sparse KKT system"
abstract type AbstractSparseKKTSystem{T, MT} <: AbstractKKTSystem{T, MT} end

#=
    Templates
=#
"Initialize KKT system with default values."
function initialize! end

"Assemble KKT matrix."
function build_kkt! end

"Compress Hessian matrix."
function compress_hessian! end

"Compress Jacobian matrix."
function compress_jacobian! end

"Get KKT system"
function get_kkt end

"Get Jacobian matrix"
function get_jacobian end

"Get Hessian matrix"
function get_hessian end

"Multiply with Jacobian"
function jtprod! end

"Regularize values in the diagonal of the KKT system."
function regularize_diagonal! end

"Set scaling of Jacobian"
function set_jacobian_scaling! end

"Return true if KKT system is reduced."
function is_reduced end

"Nonzero in Jacobian"
function nnz_jacobian end

"Nonzero in Hessian"
function nnz_kkt end

"Dense Jacobian callback"
function jac_dense! end

"Dense Hessian callback"
function hess_dense! end

#=
    Generic functions
=#
function initialize!(kkt::AbstractKKTSystem)
    fill!(kkt.pr_diag, 1.0)
    fill!(kkt.du_diag, 0.0)
    fill!(kkt.hess, 0.0)
end

function regularize_diagonal!(kkt::AbstractKKTSystem, primal, dual)
    kkt.pr_diag .+= primal
    kkt.du_diag .= .-dual
end

Base.size(kkt::AbstractKKTSystem) = size(kkt.aug_com)

# Getters
get_kkt(kkt::AbstractKKTSystem) = kkt.aug_com
get_jacobian(kkt::AbstractKKTSystem) = kkt.jac
get_hessian(kkt::AbstractKKTSystem) = kkt.hess
get_raw_jacobian(kkt::AbstractKKTSystem) = kkt.jac_raw

# Fix variable treatment
function treat_fixed_variable!(kkt::AbstractKKTSystem{T, MT}) where {T, MT<:SparseMatrixCSC{T, Int32}}
    length(kkt.ind_fixed) == 0 && return
    aug = kkt.aug_com

    fixed_aug_diag = view(aug.nzval, aug.colptr[kkt.ind_fixed])
    fixed_aug_diag .= 1.0
    fixed_aug = view(aug.nzval, kkt.ind_aug_fixed)
    fixed_aug .= 0.0
    return
end
function treat_fixed_variable!(kkt::AbstractKKTSystem{T, MT}) where {T, MT<:Matrix{T}}
    length(kkt.ind_fixed) == 0 && return
    aug = kkt.aug_com
    @inbounds for i in kkt.ind_fixed
        aug[i, :] .= 0.0
        aug[:, i] .= 0.0
        aug[i, i]  = 1.0
    end
end

#=
    SparseKKTSystem
=#

# Generic function for sparse KKT systems
function build_kkt!(kkt::AbstractSparseKKTSystem{T, MT}) where {T, MT<:Matrix{T}}
    copyto!(kkt.aug_com, kkt.aug_raw)
    treat_fixed_variable!(kkt)
end

function build_kkt!(kkt::AbstractSparseKKTSystem{T, MT}) where {T, MT<:SparseMatrixCSC{T, Int32}}
    transfer!(kkt.aug_com, kkt.aug_raw, kkt.aug_csc_map)
    treat_fixed_variable!(kkt)
end

function set_jacobian_scaling!(kkt::AbstractSparseKKTSystem{T, MT}, constraint_scaling::AbstractVector) where {T, MT}
    nnzJ = length(kkt.jac)::Int
    @inbounds for i in 1:nnzJ
        index = kkt.jac_raw.I[i]
        kkt.jacobian_scaling[i] = constraint_scaling[index]
    end
end

function compress_jacobian!(kkt::AbstractSparseKKTSystem{T, MT}) where {T, MT<:SparseMatrixCSC{T, Int32}}
    ns = length(kkt.ind_ineq)
    kkt.jac[end-ns+1:end] .= -1.0
    kkt.jac .*= kkt.jacobian_scaling # scaling
    transfer!(kkt.jac_com, kkt.jac_raw, kkt.jac_csc_map)
end

function compress_jacobian!(kkt::AbstractSparseKKTSystem{T, MT}) where {T, MT<:Matrix{T}}
    ns = length(kkt.ind_ineq)
    kkt.jac[end-ns+1:end] .= -1.0
    kkt.jac .*= kkt.jacobian_scaling # scaling
    copyto!(kkt.jac_com, kkt.jac_raw)
end

function compress_hessian!(kkt::AbstractSparseKKTSystem)
    nothing
end

function mul!(y::AbstractVector, kkt::AbstractSparseKKTSystem, x::AbstractVector)
    mul!(y, Symmetric(kkt.aug_com, :L), x)
end

function jtprod!(y::AbstractVector, kkt::AbstractSparseKKTSystem, x::AbstractVector)
    mul!(y, kkt.jac_com', x)
end

nnz_jacobian(kkt::AbstractSparseKKTSystem) = nnz(kkt.jac_raw)
nnz_kkt(kkt::AbstractSparseKKTSystem) = nnz(kkt.aug_raw)

#=
    SparseKKTSystem
=#

struct SparseKKTSystem{T, MT} <: AbstractSparseKKTSystem{T, MT}
    hess::StrideOneVector{T}
    jac::StrideOneVector{T}
    pr_diag::StrideOneVector{T}
    du_diag::StrideOneVector{T}
    # Augmented system
    aug_raw::SparseMatrixCOO{T,Int32}
    aug_com::MT
    aug_csc_map::Union{Nothing, Vector{Int}}
    # Jacobian
    jac_raw::SparseMatrixCOO{T,Int32}
    jac_com::MT
    jac_csc_map::Union{Nothing, Vector{Int}}
    # Info
    ind_ineq::Vector{Int}
    ind_fixed::Vector{Int}
    ind_aug_fixed::Vector{Int}
    jacobian_scaling::Vector{T}
end

function SparseKKTSystem{T, MT}(
    n::Int, m::Int, ind_ineq::Vector{Int}, ind_fixed::Vector{Int},
    hess_sparsity_I, hess_sparsity_J,
    jac_sparsity_I, jac_sparsity_J,
) where {T, MT}
    n_jac = length(jac_sparsity_I)
    n_hess = length(hess_sparsity_I)

    aug_vec_length = n+m
    aug_mat_length = n+m+n_hess+n_jac

    I = Vector{Int32}(undef, aug_mat_length)
    J = Vector{Int32}(undef, aug_mat_length)
    V = Vector{T}(undef, aug_mat_length)
    fill!(V, 0.0)  # Need to initiate V to avoid NaN

    offset = n+n_jac+n_hess+m

    I[1:n] .= 1:n
    I[n+1:n+n_hess] = hess_sparsity_I
    I[n+n_hess+1:n+n_hess+n_jac] .= (jac_sparsity_I.+n)
    I[n+n_hess+n_jac+1:offset] .= (n+1:n+m)

    J[1:n] .= 1:n
    J[n+1:n+n_hess] = hess_sparsity_J
    J[n+n_hess+1:n+n_hess+n_jac] .= jac_sparsity_J
    J[n+n_hess+n_jac+1:offset] .= (n+1:n+m)

    pr_diag = view(V, 1:n)
    du_diag = view(V, n_jac+n_hess+n+1:n_jac+n_hess+n+m)

    hess = view(V, n+1:n+n_hess)
    jac = view(V, n_hess+n+1:n_hess+n+n_jac)

    aug_raw = SparseMatrixCOO(aug_vec_length,aug_vec_length,I,J,V)
    jac_raw = SparseMatrixCOO(m,n,jac_sparsity_I,jac_sparsity_J,jac)

    aug_com = MT(aug_raw)
    jac_com = MT(jac_raw)

    aug_csc_map = get_mapping(aug_com, aug_raw)
    jac_csc_map = get_mapping(jac_com, jac_raw)

    ind_aug_fixed = if isa(aug_com, SparseMatrixCSC)
        _get_fixed_variable_index(aug_com, ind_fixed)
    else
        zeros(Int, 0)
    end
    jac_scaling = ones(T, n_jac)

    return SparseKKTSystem{T, MT}(
        hess, jac, pr_diag, du_diag,
        aug_raw, aug_com, aug_csc_map,
        jac_raw, jac_com, jac_csc_map,
        ind_ineq, ind_fixed, ind_aug_fixed, jac_scaling,
    )
end

# Build KKT system directly from AbstractNLPModel
function SparseKKTSystem{T, MT}(nlp::AbstractNLPModel, ind_cons=get_index_constraints(nlp)) where {T, MT}
    n_slack = length(ind_cons.ind_ineq)
    # Deduce KKT size.
    n = get_nvar(nlp) + n_slack
    m = get_ncon(nlp)
    # Evaluate sparsity pattern
    jac_I = Vector{Int32}(undef, get_nnzj(nlp))
    jac_J = Vector{Int32}(undef, get_nnzj(nlp))
    jac_structure!(nlp,jac_I, jac_J)

    hess_I = Vector{Int32}(undef, get_nnzh(nlp))
    hess_J = Vector{Int32}(undef, get_nnzh(nlp))
    hess_structure!(nlp,hess_I,hess_J)

    force_lower_triangular!(hess_I,hess_J)
    # Incorporate slack's sparsity pattern
    append!(jac_I, ind_cons.ind_ineq)
    append!(jac_J, get_nvar(nlp)+1:get_nvar(nlp)+n_slack)

    return SparseKKTSystem{T, MT}(
        n, m, ind_cons.ind_ineq, ind_cons.ind_fixed,
        hess_I, hess_J, jac_I, jac_J,
    )
end

is_reduced(::SparseKKTSystem) = true

#=
    SparseUnreducedKKTSystem
=#

struct SparseUnreducedKKTSystem{T, MT} <: AbstractSparseKKTSystem{T, MT}
    hess::StrideOneVector{T}
    jac::StrideOneVector{T}
    pr_diag::StrideOneVector{T}
    du_diag::StrideOneVector{T}

    l_diag::StrideOneVector{T}
    u_diag::StrideOneVector{T}
    l_lower::StrideOneVector{T}
    u_lower::StrideOneVector{T}

    aug_raw::SparseMatrixCOO{T,Int32}
    aug_com::MT
    aug_csc_map::Union{Nothing, Vector{Int}}

    jac_raw::SparseMatrixCOO{T,Int32}
    jac_com::MT
    jac_csc_map::Union{Nothing, Vector{Int}}
    ind_ineq::Vector{Int}
    ind_fixed::Vector{Int}
    ind_aug_fixed::Vector{Int}
    jacobian_scaling::Vector{T}
end

function SparseUnreducedKKTSystem{T, MT}(
    n::Int, m::Int, nlb::Int, nub::Int, ind_ineq, ind_fixed,
    hess_sparsity_I, hess_sparsity_J,
    jac_sparsity_I, jac_sparsity_J,
    ind_lb, ind_ub,
) where {T, MT}
    n_jac = length(jac_sparsity_I)
    n_hess = length(hess_sparsity_I)

    aug_mat_length = n + m + n_hess + n_jac + 2*nlb + 2*nub
    aug_vec_length = n+m+nlb+nub

    I = Vector{Int32}(undef, aug_mat_length)
    J = Vector{Int32}(undef, aug_mat_length)
    V = zeros(aug_mat_length)

    offset = n+n_jac+n_hess+m

    I[1:n] .= 1:n
    I[n+1:n+n_hess] = hess_sparsity_I
    I[n+n_hess+1:n+n_hess+n_jac].=(jac_sparsity_I.+n)
    I[n+n_hess+n_jac+1:offset].=(n+1:n+m)

    J[1:n] .= 1:n
    J[n+1:n+n_hess] = hess_sparsity_J
    J[n+n_hess+1:n+n_hess+n_jac].=jac_sparsity_J
    J[n+n_hess+n_jac+1:offset].=(n+1:n+m)

    I[offset+1:offset+nlb] .= (1:nlb).+(n+m)
    I[offset+nlb+1:offset+2nlb] .= (1:nlb).+(n+m)
    I[offset+2nlb+1:offset+2nlb+nub] .= (1:nub).+(n+m+nlb)
    I[offset+2nlb+nub+1:offset+2nlb+2nub] .= (1:nub).+(n+m+nlb)
    J[offset+1:offset+nlb] .= (1:nlb).+(n+m)
    J[offset+nlb+1:offset+2nlb] .= ind_lb
    J[offset+2nlb+1:offset+2nlb+nub] .= (1:nub).+(n+m+nlb)
    J[offset+2nlb+nub+1:offset+2nlb+2nub] .= ind_ub

    pr_diag = view(V,1:n)
    du_diag = view(V,n_jac+n_hess+n+1:n_jac+n_hess+n+m)

    l_diag = view(V,offset+1:offset+nlb)
    l_lower= view(V,offset+nlb+1:offset+2nlb)
    u_diag = view(V,offset+2nlb+1:offset+2nlb+nub)
    u_lower= view(V,offset+2nlb+nub+1:offset+2nlb+2nub)

    hess = view(V, n+1:n+n_hess)
    jac = view(V, n_hess+n+1:n_hess+n+n_jac)

    aug_raw = SparseMatrixCOO(aug_vec_length,aug_vec_length,I,J,V)
    jac_raw = SparseMatrixCOO(m,n,jac_sparsity_I,jac_sparsity_J,jac)

    aug_com = MT(aug_raw)
    jac_com = MT(jac_raw)

    aug_csc_map = get_mapping(aug_com, aug_raw)
    jac_csc_map = get_mapping(jac_com, jac_raw)

    jac_scaling = ones(T, n_jac)

    ind_aug_fixed = if isa(aug_com, SparseMatrixCSC)
        _get_fixed_variable_index(aug_com, ind_fixed)
    else
        zeros(Int, 0)
    end

    return SparseUnreducedKKTSystem{T, MT}(
        hess, jac, pr_diag, du_diag,
        l_diag, u_diag, l_lower, u_lower,
        aug_raw, aug_com, aug_csc_map,
        jac_raw, jac_com, jac_csc_map,
        ind_ineq, ind_fixed, ind_aug_fixed, jac_scaling,
    )
end

function SparseUnreducedKKTSystem{T, MT}(nlp::AbstractNLPModel, ind_cons=get_index_constraints(nlp)) where {T, MT}
    n_slack = length(ind_cons.ind_ineq)
    nlb = length(ind_cons.ind_lb)
    nub = length(ind_cons.ind_ub)
    # Deduce KKT size.
    n = get_nvar(nlp) + n_slack
    m = get_ncon(nlp)
    # Evaluate sparsity pattern
    jac_I = Vector{Int32}(undef, get_nnzj(nlp))
    jac_J = Vector{Int32}(undef, get_nnzj(nlp))
    jac_structure!(nlp,jac_I, jac_J)

    hess_I = Vector{Int32}(undef, get_nnzh(nlp))
    hess_J = Vector{Int32}(undef, get_nnzh(nlp))
    hess_structure!(nlp,hess_I,hess_J)

    force_lower_triangular!(hess_I,hess_J)
    # Incorporate slack's sparsity pattern
    append!(jac_I, ind_cons.ind_ineq)
    append!(jac_J, get_nvar(nlp)+1:get_nvar(nlp)+n_slack)

    return SparseUnreducedKKTSystem{T, MT}(
        n, m, nlb, nub, ind_cons.ind_ineq, ind_cons.ind_fixed,
        hess_I, hess_J, jac_I, jac_J, ind_cons.ind_lb, ind_cons.ind_ub,
    )
end

function initialize!(kkt::SparseUnreducedKKTSystem)
    kkt.pr_diag.=1
    kkt.du_diag.=0
    kkt.hess.=0
    kkt.l_lower.=0
    kkt.u_lower.=0
    kkt.l_diag.=1
    kkt.u_diag.=1
end

is_reduced(::SparseUnreducedKKTSystem) = false

#=
    DenseKKTSystem
=#

struct DenseKKTSystem{T, VT, MT} <: AbstractKKTSystem{T, MT}
    hess::MT
    jac::MT
    pr_diag::VT
    du_diag::VT
    diag_hess::VT
    # KKT system
    aug_com::MT
    # Info
    ind_ineq::Vector{Int}
    ind_fixed::Vector{Int}
    jacobian_scaling::VT
end

function DenseKKTSystem{T, VT, MT}(n, m, ind_ineq, ind_fixed) where {T, VT, MT}
    ns = length(ind_ineq)
    hess = MT(undef, n, n)
    jac = MT(undef, m, n+ns)
    pr_diag = VT(undef, n+ns)
    du_diag = VT(undef, m)
    diag_hess = VT(undef, n)

    # If the the problem is unconstrained, then KKT system is directly equal
    # to the Hessian (+ some regularization terms)
    aug_com = if (m == 0)
        hess
    else
        MT(undef, n+ns+m, n+ns+m)
    end

    jacobian_scaling = VT(undef, m)

    # Init!
    fill!(aug_com, zero(T))
    fill!(hess,    zero(T))
    fill!(jac,     zero(T))
    fill!(pr_diag, zero(T))
    fill!(du_diag, zero(T))
    fill!(diag_hess, zero(T))
    fill!(jacobian_scaling, one(T))

    return DenseKKTSystem{T, VT, MT}(
        hess, jac, pr_diag, du_diag, diag_hess, aug_com, ind_ineq, ind_fixed, jacobian_scaling,
    )
end

function DenseKKTSystem{T, VT, MT}(nlp::AbstractNLPModel, info_constraints=get_index_constraints(nlp)) where {T, VT, MT}
    return DenseKKTSystem{T, VT, MT}(
        get_nvar(nlp), get_ncon(nlp), info_constraints.ind_ineq, info_constraints.ind_fixed,
    )
end

is_reduced(::DenseKKTSystem) = true

# Special getters for Jacobian
function get_jacobian(kkt::DenseKKTSystem)
    n = size(kkt.hess, 1)
    ns = length(kkt.ind_ineq)
    return view(kkt.jac, :, 1:n)
end
get_raw_jacobian(kkt::DenseKKTSystem) = kkt.jac

nnz_jacobian(kkt::DenseKKTSystem) = length(kkt.jac)
nnz_kkt(kkt::DenseKKTSystem) = length(kkt.aug_com)

function _update_diagonal!(dest::AbstractMatrix, d1::AbstractVector, d2::AbstractVector)
    n = length(d1)
    for i in 1:n
        dest[i, i] = d1[i] + d2[i]
    end
end

function _build_dense_kkt_system!(dest, hess, jac, pr_diag, du_diag, diag_hess, n, m, ns)
    # Transfer Hessian
    for i in 1:n, j in 1:i
        if i == j
            dest[i, i] = pr_diag[i] + diag_hess[i]
        else
            dest[i, j] = hess[i, j]
            dest[j, i] = hess[j, i]
        end
    end
    # Transfer slack diagonal
    for i in 1:ns
        dest[i+n, i+n] = pr_diag[i+n]
    end
    # Transfer Jacobian
    for i in 1:m, j in 1:(n+ns)
        dest[i + n + ns, j] = jac[i, j]
        dest[j, i + n + ns] = jac[i, j]
    end
    # Transfer dual regularization
    for i in 1:m
        dest[i + n + ns, i + n + ns] = du_diag[i]
    end
end

function build_kkt!(kkt::DenseKKTSystem{T, VT, MT}) where {T, VT, MT}
    n = size(kkt.hess, 1)
    m = size(kkt.jac, 1)
    ns = length(kkt.ind_ineq)
    if m == 0 # If problem is unconstrained, just need to update the diagonal
        _update_diagonal!(kkt.aug_com, kkt.diag_hess, kkt.pr_diag)
    else # otherwise, we update the full matrix
        _build_dense_kkt_system!(kkt.aug_com, kkt.hess, kkt.jac, kkt.pr_diag, kkt.du_diag, kkt.diag_hess, n, m, ns)
    end
    treat_fixed_variable!(kkt)
end

function compress_jacobian!(kkt::DenseKKTSystem{T, VT, MT}) where {T, VT, MT}
    m = size(kkt.jac, 1)
    n = size(kkt.hess, 1)
    # Add slack indexes
    for i in kkt.ind_ineq
        kkt.jac[i, i+n] = -one(T)
    end
    # Scale
    kkt.jac .*= kkt.jacobian_scaling
    return
end

function compress_hessian!(kkt::DenseKKTSystem)
    # Transfer diagonal term for future regularization
    diag!(kkt.diag_hess, kkt.hess)
end

function mul!(y::AbstractVector, kkt::DenseKKTSystem, x::AbstractVector)
    mul!(y, kkt.aug_com, x)
end

function jtprod!(y::AbstractVector, kkt::DenseKKTSystem, x::AbstractVector)
    mul!(y, kkt.jac', x)
end

function set_jacobian_scaling!(kkt::DenseKKTSystem, constraint_scaling::AbstractVector)
    copyto!(kkt.jacobian_scaling, constraint_scaling)
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_checklist

"""
    dcc_checklist(;kwargs...)

A Checklist component.
Checklist is a component that encapsulates several checkboxes.
The values and labels of the checklist are specified in the `options`
property and the checked items are specified with the `value` property.
Each checkbox is rendered as an input with a surrounding label.
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `options` (optional): An array of options. options has the following type: Array of lists containing elements 'label', 'value', 'disabled'.
Those elements have the following types:
  - `label` (String | Real; required): The checkbox's label
  - `value` (String | Real; required): The value of the checkbox. This value
corresponds to the items specified in the
`value` property.
  - `disabled` (Bool; optional): If true, this checkbox is disabled and can't be clicked on.s
- `value` (Array of String | Reals; optional): The currently selected value
- `className` (String; optional): The class of the container (div)
- `style` (Dict; optional): The style of the container (div)
- `inputStyle` (Dict; optional): The style of the <input> checkbox element
- `inputClassName` (String; optional): The class of the <input> checkbox element
- `labelStyle` (Dict; optional): The style of the <label> that wraps the checkbox input
 and the option's label
- `labelClassName` (String; optional): The class of the <label> that wraps the checkbox input
 and the option's label
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `persistence` (Bool | String | Real; optional): Used to allow user interactions in this component to be persisted when
the component - or the page - is refreshed. If `persisted` is truthy and
hasn't changed from its previous value, a `value` that the user has
changed while using the app will keep that change, as long as
the new `value` also matches what was given originally.
Used in conjunction with `persistence_type`.
- `persisted_props` (Array of a value equal to: 'value's; optional): Properties whose user interactions will persist after refreshing the
component or the page. Since only `value` is allowed this prop can
normally be ignored.
- `persistence_type` (a value equal to: 'local', 'session', 'memory'; optional): Where persisted user changes will be stored:
memory: only kept in memory, reset on page refresh.
local: window.localStorage, data is kept after the browser quit.
session: window.sessionStorage, data is cleared once the browser quit.
"""
function dcc_checklist(; kwargs...)
        available_props = Symbol[:id, :options, :value, :className, :style, :inputStyle, :inputClassName, :labelStyle, :labelClassName, :loading_state, :persistence, :persisted_props, :persistence_type]
        wild_props = Symbol[]
        return Component("dcc_checklist", "Checklist", "dash_core_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
function f()
    n = 1
    if n == 1
        print("a")
    end
end

n = 0
f()
if n == 0
    print("b")
end
print("\n")

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function rangesearch!{T}(metric::Function, neighbors::NeighborList{T}, node::Node{T}, point_index::Int, r::T)
    node.index == 0 && return

    d = metric(node.index, point_index)::T
    d <= r && add_neighbor!(neighbors, node.index, d)

    node.isleaf && return

    μ = node.radius
    if d + μ <= r
        add_neighbors!(metric, neighbors, node.inside, point_index)
    elseif d - μ <= r
        rangesearch!(metric, neighbors, node.inside, point_index, r)
    end

    if μ - d <= r
        rangesearch!(metric, neighbors, node.outside, point_index, r)
    end
end

function rangesearch_sq!{T}(metric::Function, neighbors::NeighborList{T}, node::Node{T}, point_index::Int, r::T)
    node.index == 0 && return

    d = metric(node.index, point_index)::T
    d <= r && add_neighbor!(neighbors, node.index, d)

    node.isleaf && return

    μ = node.radius
    # d + μ <= r
    if !sumsqrtgt(d, μ, r)
        add_neighbors!(metric, neighbors, node.inside, point_index)
    # d - μ <= r
    elseif sumsqrtgt(r, μ, d)
        rangesearch_sq!(metric, neighbors, node.inside, point_index, r)
    end

    # μ - d <= r
    if sumsqrtgt(r, d, μ)
        rangesearch_sq!(metric, neighbors, node.outside, point_index, r)
    end
end

# Searches using a possible *new* metric
function rangesearch!{T}(metric::Function, neighbors::NeighborList{T}, tree::VPTree, point_index::Int, r::T)
    clear!(neighbors)
    rangesearch!(metric, neighbors, tree.root, point_index, r)
    return neighbors
end

function rangesearch!{T}(neighbors::NeighborList{T}, tree::VPTree, point_index::Int, r::T)
    clear!(neighbors)
    rangesearch!(tree.metric, neighbors, tree.root, point_index, r)
    return neighbors
end

function rangesearch_sq!{T}(neighbors::NeighborList{T}, tree::VPTree, point_index::Int, r::T)
    clear!(neighbors)
    rangesearch_sq!(tree.metric, neighbors, tree.root, point_index, r)
    return neighbors
end

function rangesearch{T}(tree::VPTree, point_index::Int, r::T)
    neighbors = NeighborList{T}(tree.num_points)
    rangesearch!(neighbors, tree, point_index, r)
end

function rangesearch_sq{T}(tree::VPTree, point_index::Int, r::T)
    neighbors = NeighborList{T}(tree.num_points)
    rangesearch_sq!(neighbors, tree, point_index, r)
end


text:
close("all")
transD_GP.getchi2forall(opt)
ax = gcf().axes;
χ² = length(F.d_log10_ρ)
ax[2].plot(xlim(), [χ², χ²], "--", color="gray")
ax[2].set_ylim(χ²-10, χ²+10)
##
opt.xall[:] .= zall
transD_GP.plot_posterior(F, opt, burninfrac=0.5, figsize=(10,6), qp1=0.05, qp2=0.95, nbins=50, vmaxpc=1.0)
ax = gcf().axes
ax[1].step(log10.(ρ[2:end]), z[2:end], color="w", linewidth=2)
ax[1].step(log10.(ρ[2:end]), z[2:end], color="g", linestyle="--")

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
#!/usr/bin/env julia

# This is based on the runtests.jl file of MPI.jl.

using MPI: MPI, mpiexec
using InteractiveUtils: versioninfo

# Load test packages to trigger precompilation
using PencilFFTs

test_files = [
    "taylor_green.jl",
    "rfft.jl",
    "transforms.jl",
]

# Make sure that example files run.
example_dir = joinpath("..", "examples")
example_files = joinpath.(
    example_dir,
    filter(fname -> splitext(fname)[2] == ".jl", readdir(example_dir))
)

Nproc = let N = get(ENV, "JULIA_MPI_TEST_NPROCS", nothing)
    N === nothing ? clamp(Sys.CPU_THREADS, 4, 6) : parse(Int, N)
end
files = [test_files..., example_files...]

println()
versioninfo()
println("\n", MPI.MPI_LIBRARY_VERSION_STRING, "\n")

for fname in files
    @info "Running $fname with $Nproc processes..."
    mpiexec() do cmd
        # Disable precompilation to prevent race conditions when loading
        # packages.
        run(`$cmd -n $Nproc $(Base.julia_cmd()) $fname`)
    end
    println()
end


text:
module Example

export MyTypeA, MyTypeB
export greet

abstract type AbstractMyType end

struct MyTypeA <: AbstractMyType end

struct MyTypeB <: AbstractMyType end

greet(::Type{MyTypeA}) = "Hello, this is MyTypeA."
greet(::Type{MyTypeB}) = "Hi, this is MyTypeB!"

end # module


text:
using Transformers
using Transformers.Basic
using Test

const tests = [
    "transformer",
    "embed",
]

@test isempty(detect_ambiguities(Transformers, Transformers.GenerativePreTrain, Transformers.Basic, Transformers.Datasets))

@testset "Transformers" begin
    for t in tests
        fp = joinpath(dirname(@__FILE__), "test_$t.jl")
        println("$fp ...")
        include(fp)
    end
end


text:
# APT module relying on sudo aptitude commands
module APT

export Package, readpkgs, writepkgs, rmpkgs!, addpkgs!
export update, install, remove, clean

immutable Package
    name::ASCIIString
end

pkglt(x, y) = x.name < y.name

sort!(pkgs::Vector{Package}) = Base.sort!(pkgs, lt=pkglt)

## Package file management

function readpkg(stream::IOStream)
    line = chomp(readline(stream))

    if line[1] == '+'
        (:install, Package(line[2:end]))
    elseif line[1] == '-'
        (:remove, Package(line[2:end]))
    else
        error("Bad package line.")
    end
end

function readpkgs(stream::IOStream,
                  installpkgs::Vector{Package},
                  removepkgs::Vector{Package})
    while !eof(stream)
        op, pkg = readpkg(stream)

        if op == :install
            push!(installpkgs, pkg)
        elseif op == :remove
            push!(removepkgs, pkg)
        else
            error("Bad package op.")
        end
    end

    (sort!(installpkgs), sort!(removepkgs))
end

function readpkgs(filepath::ASCIIString)
    installpkgs = Package[]
    removepkgs = Package[]

    try
        file = open(filepath)
        readpkgs(file, installpkgs, removepkgs)
        close(file)
        (installpkgs, removepkgs)
    catch
        error("Bad file path.")
    end
end

function writepkg(stream::IOStream, op::Symbol, pkg::Package)
    if op == :install
        write(stream, "+$(pkg.name)\n")
    elseif op == :remove
        write(stream, "-$(pkg.name)\n")
    else
        error("Bad package op.")
    end
end

function writepkgs(stream::IOStream,
                   installpkgs::Vector{Package},
                   removepkgs::Vector{Package})
    map(pkg -> writepkg(stream, :install, pkg), sort!(installpkgs))
    map(pkg -> writepkg(stream, :remove, pkg), sort!(removepkgs))
end

function writepkgs(filepath::ASCIIString,
                   installpkgs::Vector{Package},
                   removepkgs::Vector{Package})
        file = open(filepath, "w")
        writepkgs(file, installpkgs, removepkgs)
        close(file)
end

rmpkg!(pkg::Package, frompkgs::Vector{Package}) =
    filter!(p -> p.name != pkg.name, frompkgs)

rmpkgs!(pkgs::Vector{Package}, frompkgs::Vector{Package}) =
    map(pkg -> rmpkg!(pkg, frompkgs), pkgs)

rmpkgs!(pkgnames::Vector{ASCIIString}, frompkgs::Vector{Package}) =
    map(pkgname -> rmpkg!(Package(pkgname), frompkgs), pkgnames)

function in(pkg::Package, pkgs::Vector{Package})
    if length(pkgs) == 0
        false
    elseif length(pkgs) == 1
        pkg.name == pkgs[1].name
    else
        in(pkg, pkgs[1:1]) || in(pkg, pkgs[2:end])
    end
end

function addpkg!(pkg::Package,
                 topkgs::Vector{Package},
                 frompkgs::Vector{Package})
    if !in(pkg, topkgs)
        push!(topkgs, pkg)
        rmpkg!(pkg, frompkgs)
    end
end

addpkgs!(pkgs::Vector{Package},
         topkgs::Vector{Package},
         frompkgs::Vector{Package}) =
    map(pkg -> addpkg!(pkg, topkgs, frompkgs), pkgs)

addpkgs!(pkgnames::Vector{ASCIIString},
         topkgs::Vector{Package},
         frompkgs::Vector{Package}) =
    map(pkgname -> addpkg!(Package(pkgname), topkgs, frompkgs), pkgnames)

## aptitude commands

update() = run(`sudo aptitude update`)

install(pkgs::Vector{Package}) =
    run(`sudo aptitude install $([pkg.name for pkg in pkgs])`)
install(pkgnames::Vector{ASCIIString}) = install(map(Package, pkgnames))

remove(pkgs::Vector{Package}) =
    run(`sudo aptitude --purge remove $([pkg.name for pkg in pkgs])`)
remove(pkgnames::Vector{ASCIIString}) = remove(map(Package, pkgnames))

upgrade() = run(`sudo aptitude upgrade`)

clean() = run(`sudo aptitude clean`)

end # module


text:
## Imports

using Aqua
using ImplicitDifferentiation
using JET
using Random
using Test

## Test sets

@testset verbose = true "ImplicitDifferentiation.jl" begin
    @testset verbose = true "Code quality" begin
        @testset verbose = true "JET" begin
            jet_report = JET.report_package(ImplicitDifferentiation)
            @test string(jet_report) == "No errors detected\n"
        end
        @testset verbose = true "Aqua" begin
            Aqua.test_all(ImplicitDifferentiation)
        end
    end
    @testset verbose = true "Unconstrained optimization" begin
        include("1_unconstrained_optimization.jl")
    end
    @testset verbose = true "Constrained optimization" begin
        include("2_constrained_optimization.jl")
    end
    @testset verbose = true "Entropy regularized optimal transport" begin
        include("3_optimal_transport.jl")
    end
    @testset verbose = true "Custom structs" begin
        include("4_struct.jl")
    end
end


text:
"""
$(TYPEDEF)

Store both grid and density for KDE over the real line.

Reading the fields directly is part of the API, and

```julia
sum(density) * step(x) * step(y) ≈ 1
```

# Fields

$(FIELDS)
"""
mutable struct BivariateKDE{Rx<:AbstractRange,Ry<:AbstractRange} <: AbstractKDE
    "First coordinate of gridpoints for evaluating the density."
    x::Rx
    "Second coordinate of gridpoints for evaluating the density."
    y::Ry
    "Kernel density at corresponding gridpoints `Tuple.(x, permutedims(y))`."
    density::Matrix{Float64}
end

function kernel_dist(::Type{D},w::Tuple{Real,Real}) where D<:UnivariateDistribution
    kernel_dist(D,w[1]), kernel_dist(D,w[2])
end
function kernel_dist(::Type{Tuple{Dx, Dy}},w::Tuple{Real,Real}) where {Dx<:UnivariateDistribution,Dy<:UnivariateDistribution}
    kernel_dist(Dx,w[1]), kernel_dist(Dy,w[2])
end

const DataTypeOrUnionAll = Union{DataType, UnionAll}

# this function provided for backwards compatibility, though it doesn't have the type restrictions
# to ensure that the given tuple only contains univariate distributions
function kernel_dist(d::Tuple{DataTypeOrUnionAll, DataTypeOrUnionAll}, w::Tuple{Real,Real})
    kernel_dist(d[1],w[1]), kernel_dist(d[2],w[2])
end

# TODO: there are probably better choices.
function default_bandwidth(data::Tuple{RealVector,RealVector})
    default_bandwidth(data[1]), default_bandwidth(data[2])
end

# tabulate data for kde
function tabulate(data::Tuple{RealVector, RealVector}, midpoints::Tuple{Rx, Ry},
        weights::Weights = default_weights(data)) where {Rx<:AbstractRange,Ry<:AbstractRange}
    xdata, ydata = data
    ndata = length(xdata)
    length(ydata) == ndata || error("data vectors must be of same length")

    xmid, ymid = midpoints
    nx, ny = length(xmid), length(ymid)
    sx, sy = step(xmid), step(ymid)

    # Set up a grid for discretized data
    grid = zeros(Float64, nx, ny)
    ainc = 1.0 / (sum(weights)*(sx*sy)^2)

    # weighted discretization (cf. Jones and Lotwick)
    for i in 1:length(xdata)
        x = xdata[i]
        y = ydata[i]
        kx, ky = searchsortedfirst(xmid,x), searchsortedfirst(ymid,y)
        jx, jy = kx-1, ky-1
        if 1 <= jx <= nx-1 && 1 <= jy <= ny-1
            grid[jx,jy] += (xmid[kx]-x)*(ymid[ky]-y)*ainc*weights[i]
            grid[kx,jy] += (x-xmid[jx])*(ymid[ky]-y)*ainc*weights[i]
            grid[jx,ky] += (xmid[kx]-x)*(y-ymid[jy])*ainc*weights[i]
            grid[kx,ky] += (x-xmid[jx])*(y-ymid[jy])*ainc*weights[i]
        end
    end

    # returns an un-convolved KDE
    BivariateKDE(xmid, ymid, grid)
end

# convolution with product distribution of two univariates distributions
function conv(k::BivariateKDE, dist::Tuple{UnivariateDistribution,UnivariateDistribution})
    # Transform to Fourier basis
    Kx, Ky = size(k.density)
    ft = rfft(k.density)

    distx, disty = dist

    # Convolve fft with characteristic function of kernel
    cx = -twoπ/(step(k.x)*Kx)
    cy = -twoπ/(step(k.y)*Ky)
    for j = 0:size(ft,2)-1
        for i = 0:size(ft,1)-1
            ft[i+1,j+1] *= cf(distx,i*cx)*cf(disty,min(j,Ky-j)*cy)
        end
    end
    dens = irfft(ft, Kx)

    for i = 1:length(dens)
        dens[i] = max(0.0,dens[i])
    end

    # Invert the Fourier transform to get the KDE
    BivariateKDE(k.x, k.y, dens)
end

const BivariateDistribution = Union{MultivariateDistribution,Tuple{UnivariateDistribution,UnivariateDistribution}}

default_weights(data::Tuple{RealVector, RealVector}) = UniformWeights(length(data[1]))

function kde(data::Tuple{RealVector, RealVector}, weights::Weights, midpoints::Tuple{Rx, Ry},
        dist::BivariateDistribution) where {Rx<:AbstractRange,Ry<:AbstractRange}
    k = tabulate(data, midpoints, weights)
    conv(k,dist)
end

function kde(data::Tuple{RealVector, RealVector}, dist::BivariateDistribution;
             boundary::Tuple{Tuple{Real,Real}, Tuple{Real,Real}} = (kde_boundary(data[1],std(dist[1])),
                                                     kde_boundary(data[2],std(dist[2]))),
             npoints::Tuple{Int,Int}=(256,256),
             weights::Weights = default_weights(data))

    xmid = kde_range(boundary[1],npoints[1])
    ymid = kde_range(boundary[2],npoints[2])

    kde(data,weights,(xmid,ymid),dist)
end

function kde(data::Tuple{RealVector, RealVector}, midpoints::Tuple{Rx, Ry};
             bandwidth=default_bandwidth(data), kernel=Normal,
             weights::Weights = default_weights(data)) where {Rx<:AbstractRange,Ry<:AbstractRange}

    dist = kernel_dist(kernel,bandwidth)
    kde(data,weights,midpoints,dist)
end

function kde(data::Tuple{RealVector, RealVector};
             bandwidth=default_bandwidth(data),
             kernel=Normal,
             boundary::Tuple{Tuple{Real,Real}, Tuple{Real,Real}} = (kde_boundary(data[1],bandwidth[1]),
                                                     kde_boundary(data[2],bandwidth[2])),
             npoints::Tuple{Int,Int}=(256,256),
             weights::Weights = default_weights(data))

    dist = kernel_dist(kernel,bandwidth)
    xmid = kde_range(boundary[1],npoints[1])
    ymid = kde_range(boundary[2],npoints[2])

    kde(data,weights,(xmid,ymid),dist)
end

# matrix data
function kde(data::RealMatrix,args...;kwargs...)
    size(data,2) == 2 || error("Can only construct KDE from matrices with 2 columns.")
    kde((data[:,1],data[:,2]),args...;kwargs...)
end


text:
abstract type ConstrainedOptimizer{T}  <: AbstractConstrainedOptimizer end
abstract type IPOptimizer{T} <: ConstrainedOptimizer{T} end # interior point methods


text:
# Powermethod using Tpetra
# Adapted from https://trilinos.org/docs/dev/packages/tpetra/doc/html/Tpetra_Lesson03.html

using MPI
using Trilinos

function powermethod(A, niters::Number, tolerance::Number)
  q = Tpetra.Vector(Tpetra.getDomainMap(A))
  z = Tpetra.Vector(Tpetra.getRangeMap(A))
  resid = Tpetra.Vector(Tpetra.getRangeMap(A))

  Tpetra.randomize(z)

  lambda = 0.0
  normz = 0.0
  residual = 0.0

  report_frequency = 10

  for iter in 1:niters
    normz = Tpetra.norm2(z)         # Compute the 2-norm of z
    Tpetra.scale(q, 1.0 / normz, z) # q := z / normz
    Tpetra.apply(A, q, z)           # z := A * q
    lambda = dot(q,z)               # Approx. max eigenvalue
    # Compute and report the residual norm every report_frequency
    # iterations, or if we've reached the maximum iteration count.
    if iter % report_frequency == 0 || iter == 1
      Tpetra.update(resid, 1.0, z, -lambda, q, 0.0) # z := A*q - lambda*q
      residual = Tpetra.norm2(resid)                     # 2-norm of the residual vector

      println("Iteration $iter:")
      println("- lambda = $lambda")
      println("- ||A*q - lambda*q||_2 = $residual")
    end
    if residual < tolerance
      println("Converged after $iter iterations")
      break
    elseif (iter == niters)
      println("Failed to converge after $niters iterations")
      break
    end
  end

  return lambda
end

# MPI setup
MPI.Init()
comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))

# Communication map
const num_global_indices = 50
const indexbase = 0
map = Tpetra.Map(num_global_indices, indexbase, comm)

# Number of local elements
my_rank = Teuchos.getRank(comm)
if my_rank != 0
  redirect_stdout(open("/dev/null", "w"))
  redirect_stderr(open("/dev/null", "w"))
end

num_my_elements = Int(Tpetra.getNodeNumElements(map)) # UInt otherwise
println("Number of elements for rank $my_rank is $num_my_elements")

# Construct the matrix
A = Tpetra.CrsMatrix(map, 0)
for local_row in 0:num_my_elements-1
  global_row = Tpetra.getGlobalElement(map, local_row)
  if global_row == 0
    Tpetra.insertGlobalValues(A, global_row, [global_row, global_row+1], [2.0, -1.0])
  elseif global_row == num_my_elements-1
    Tpetra.insertGlobalValues(A, global_row, [global_row-1, global_row], [-1.0, 2.0])
  else
    Tpetra.insertGlobalValues(A, global_row, [global_row-1, global_row, global_row+1], [-1.0, 2.0, -1.0])
  end
end

Tpetra.fillComplete(A)

# Number of iterations
niters = 500
# Target tolerance
tolerance = 1e-2

lambda = powermethod(A, niters, tolerance)
println("\nMaximal estimated eigenvalue: $lambda\n")

println("Increasing magnitude of A(0,0), solving again")

Tpetra.resumeFill(A)

if Tpetra.isNodeGlobalElement(Tpetra.getRowMap(A), 0)
  first_row_id = 0
  num_entries_in_row = Tpetra.getNumEntriesInGlobalRow(A, first_row_id)

  row_values = Vector{Float64}(num_entries_in_row)
  row_indices = Vector{Int}(num_entries_in_row)

  num_entries_returned = Ref(UInt(0))

  # The last 3 parameters here are output parameters. The ArrayView wraps the provided Julia array, which must be at least as large as the number of entries stored
  # num_entries_returned contains the actual number of items that is read
  Tpetra.getGlobalRowCopy(A, first_row_id, Teuchos.ArrayView(row_indices), Teuchos.ArrayView(row_values), num_entries_returned)

  for i in 1:num_entries_returned[]
    if row_indices[i] == first_row_id
      row_values[i] *= 10.0
    end
  end

  # Conversion to ArrayView is implicit here, because the arrays are read-only
  Tpetra.replaceGlobalValues(A, first_row_id, row_indices, row_values)
end

Tpetra.fillComplete(A)

lambda = powermethod(A, niters, tolerance)
println("\nMaximal estimated eigenvalue: $lambda\n")

MPI.Finalize()


text:
const monadic_nonlinear_special = [airyai, airyaiprime, airybi, airybiprime, besselj0, besselj1, bessely0, bessely1, dawson, digamma, erf, erfc, erfcinv, erfi, erfinv, erfcx, gamma, invdigamma, lgamma, trigamma]

const diadic_linear_false_false_false_special = [lbeta, beta]
# const diadic_linear_true_false_true_special = [besselk, hankelh2, bessely, besselj, besseli, polygamma, hankelh1]

for f in monadic_nonlinear_special
    @eval begin
        haslinearity(::typeof($f), ::Val{1}) = true
        linearity(::typeof($f), ::Val{1}) = Val(false)
    end
end

# linearity of a 2-arg function is:
# Val((linear11, linear22, linear12))
#
# linearIJ refers to the zeroness of d^2/dxIxJ
for f in diadic_linear_false_false_false_special
    @eval begin
        haslinearity(::typeof($f), ::Val{2}) = true
        linearity(::typeof($f), ::Val{2}) = Val((false,false,false))
    end
end

# for f in diadic_linear_true_false_true_special
#     @eval begin
#         haslinearity(::typeof($f), ::Val{2}) = true
#         linearity(::typeof($f), ::Val{2}) = Val((true,false,true))
#     end
# end


text:
type PottsModel{T} <: ClassicalSpinModel{T}
	spininfo::SpinInfo{T}
	factorW::Array{T,2}

	function PottsModel{T}(spininfo::SpinInfo{T})
		this = new{T}()
		this.spininfo = checkSpininfoForPotts(spininfo)
		initialize!(this)
		return this
	end
end

PottsModel{T}(spininfo::SpinInfo{T}) = PottsModel{T}(spininfo)
PottsModel() = testPottsModel()

#---
# functions for constructor
function checkSpininfoForPotts{T}(spininfo::SpinInfo{T})
	if spininfo.modelname != "potts"
		println("warning:spininfo.modelname does not match to the type.
		  the modelname is changed to 'potts'")
		spininfo.modelname = "potts"
	end
	return spininfo
end

function initialize!{T}(pottsmodel::PottsModel{T})
	if isSymmetricFactorization(pottsmodel)
		factorW = initializePottsFactorWsym(pottsmodel::PottsModel)
	else
		factorW = initializeFactorWasym(pottsmodel::PottsModel)
	end
	pottsmodel.factorW = factorW
end

function initializePottsFactorWsym{T}(pottsmodel::PottsModel{T})
	temp, externalfield = getEnvParameters(pottsmodel)
	numberOfState = getStates(pottsmodel)
	sigmah = 1
	revTemp = one(T)/temp	# this is K
	etoK = exp(revTemp)		# e to the K
	eKbar = etoK + sqrt( (etoK + numberOfState - one(T)) * (etoK - one(T)) )
	eGamma = getExpGamma(revTemp, externalfield)
	denominatorTerm = sqrt( numberOfState -2.0 + 2.0 * eKbar )
	factorW = ones(T, numberOfState, numberOfState)
	for sigma = 1:numberOfState, s = 1:numberOfState
		if sigma == s
			factorW[sigma, s] *= eKbar
		end
		if sigma == sigmah
			factorW[sigma, s] *= eGamma
		end
	end
	factorW = factorW ./ denominatorTerm
	return factorW
end

function getHamiltonian{T}(pottsmodel::PottsModel{T})
	externalfield = getExternalfield(pottsmodel)
	numberOfState = getStates(pottsmodel)
	sigmah = 1
	hamiltonian = zeros(T, numberOfState, numberOfState)
	for sigmai = 1:numberOfState, sigmaj = 1:numberOfState
		addterm = zero(T)
		if sigmai == sigmah
			addterm += externalfield / 4.0
		end
		if sigmaj == sigmah
			addterm += externalfield / 4.0
		end
		if sigmai == sigmaj
			addterm += one(T)
		end
		hamiltonian[sigmai, sigmaj] += addterm
	end
	return -hamiltonian
end

function getMeasureOperator{T}(spinmodel::PottsModel{T})
	measurementOperator = zeros(getStates(spinmodel))
	measurementOperator[1] = one(T)
	return measurementOperator
end

#---
# functions for PottsModel
function isPottsModel(pottsmodel::PottsModel)
	return true
end

#---
# test to make a example SpinModel
function testPottsModel(numberOfState::Int = 2)
	modelcode = "classical_potts_" * dec(numberOfState)
	factorization = "sym"
	externalfield = 1.0e-13
	temperature = 1.5
	spininfo = SpinInfo(modelcode, factorization, externalfield, temperature)
	testModel = buildSpinSystem(spininfo)
	return testModel
end


text:
@doc doc"""
    Euclidean{T<:Tuple} <: Manifold

Euclidean vector space $\mathbb R^n$.

# Constructor

    Euclidean(n)

Generate the $n$-dimensional vector space $\mathbb R^n$.

    Euclidean(n₁,n₂,...,nᵢ; field=ℝ)

Generate the vector space of $k=n_1n_2\cdot\ldots n_i$ values, i.e. the
$\mathbb F^{n_1, n_2, \ldots, n_d}$ whose
elements are interpreted as $n_1 \times,n_2\times\cdots\times n_i$ arrays.
For $d=2$ we obtain a matrix space.
The default `field=ℝ` can also be set to `field=ℂ`.
The dimension of this space is $k \dim_ℝ 𝔽$, where $\dim_ℝ 𝔽$ is the
[`real_dimension`](@ref) of the field $𝔽$.
"""
struct Euclidean{N<:Tuple,T} <: Manifold where {N,T<:AbstractNumbers} end

function Euclidean(n::Vararg{Int,N}; field::AbstractNumbers = ℝ) where {N}
    return Euclidean{Tuple{n...},field}()
end

"""
    EuclideanMetric <: RiemannianMetric

A general type for any manifold that employs the Euclidean Metric, for example
the [`Euclidean`](@ref) manifold itself, or the [`Sphere`](@ref), where every
tangent space (as a plane in the embedding) uses this metric (in the embedding).

Since the metric is independent of the field type, this metric is also used for
the Hermitian metrics, i.e. metrics that are analogous to the `EuclideanMetric`
but where the field type of the manifold is `ℂ`.

This metric is the default metric for example for the [`Euclidean`](@ref) manifold.
"""
struct EuclideanMetric <: RiemannianMetric end

det_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = one(eltype(x))

"""
    distance(M::Euclidean, x, y)

Compute the Euclidean distance between two points on the [`Euclidean`](@ref)
manifold `M`, i.e. for vectors it's just the norm of the difference, for matrices
and higher order arrays, the matrix and ternsor Frobenius norm, respectively.
"""
distance(::Euclidean, x, y) = norm(x .- y)

@doc doc"""
    exp(M::Euclidean, x, v)

Compute the exponential map on the [`Euclidean`](@ref) manifold `M` from `x` in direction
`v`, which in this case is just
````math
\exp_x v = x + v.
````
"""
exp(::Euclidean, ::Any...)

exp!(M::Euclidean, y, x, v) = (y .= x .+ v)

"""
    flat(M::Euclidean, x, w)

Transform a tangent vector into a cotangent. Since they can directly be identified in the
[`Euclidean`](@ref) case, this yields just the identity for a tangent vector `w` in the
tangent space of `x` on `M`. The result is returned also in place in `v`.
"""
flat(::Euclidean, ::Any...)

flat!(M::Euclidean, v::CoTFVector, x, w::TFVector) = copyto!(v, w)

function get_basis(M::Euclidean{<:Tuple,ℝ}, x, B::ArbitraryOrthonormalBasis)
    vecs = [_euclidean_basis_vector(x, i) for i in eachindex(x)]
    return PrecomputedOrthonormalBasis(vecs)
end
function get_basis(M::Euclidean{<:Tuple,ℂ}, x, B::ArbitraryOrthonormalBasis)
    vecs = [_euclidean_basis_vector(x, i) for i in eachindex(x)]
    return PrecomputedOrthonormalBasis([vecs; im * vecs])
end
function get_basis(M::Euclidean, x, B::DiagonalizingOrthonormalBasis)
    vecs = get_basis(M, x, ArbitraryOrthonormalBasis()).vectors
    kappas = zeros(real(eltype(x)), manifold_dimension(M))
    return PrecomputedDiagonalizingOrthonormalBasis(vecs, kappas)
end

function get_coordinates(M::Euclidean{<:Tuple,ℝ}, x, v, B::ArbitraryOrDiagonalizingBasis)
    S = representation_size(M)
    PS = prod(S)
    return reshape(v, PS)
end
function get_coordinates(M::Euclidean{<:Tuple,ℂ}, x, v, B::ArbitraryOrDiagonalizingBasis)
    S = representation_size(M)
    PS = prod(S)
    return [reshape(real(v), PS); reshape(imag(v), PS)]
end

function get_vector(M::Euclidean{<:Tuple,ℝ}, x, v, B::ArbitraryOrDiagonalizingBasis)
    S = representation_size(M)
    return reshape(v, S)
end
function get_vector(M::Euclidean{<:Tuple,ℂ}, x, v, B::ArbitraryOrDiagonalizingBasis)
    S = representation_size(M)
    N = div(length(v), 2)
    return reshape(v[1:N] + im * v[N+1:end], S)
end

function hat(M::Euclidean{N,ℝ}, x, vⁱ) where {N}
    return reshape(vⁱ, representation_size(TangentBundleFibers(M)))
end

hat!(::Euclidean{N,ℝ}, v, x, vⁱ) where {N} = copyto!(v, vⁱ)

@doc doc"""
    injectivity_radius(M::Euclidean)

Return the injectivity radius on the [`Euclidean`](@ref) `M`, which is $\infty$.
"""
injectivity_radius(::Euclidean) = Inf

@doc doc"""
    inner(M::Euclidean, x, v, w)

Compute the inner product on the [`Euclidean`](@ref) `M`, which is just
the inner product on the real-valued or complex valued vector space
of arrays (or tensors) of size $n_1\times n_2 \times \cdots \times n_i$, i.e.

````math
g_x(v,w) = \sum_{k\in I} \overline{v}_{k} w_{k},
````
where $I$ is the set of integer vectors $k\in\mathbb N^i$, such that for all
$1 \leq j \leq i$ it holds $1\leq k_j \leq n_j$.

For the special case of $i\leq 2$, i.e. matrices and vectors, this simplifies to
````math
g_x(v,w) = w^{\mathrm{H}}v,
````
where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
"""
inner(::Euclidean, ::Any...)
@inline inner(::Euclidean, x, v, w) = dot(v, w)
@inline inner(::MetricManifold{<:Manifold,EuclideanMetric}, x, v, w) = dot(v, w)

inverse_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = local_metric(M, x)

is_default_metric(::Euclidean, ::EuclideanMetric) = Val(true)

function local_metric(::MetricManifold{<:Manifold,EuclideanMetric}, x)
    return Diagonal(ones(SVector{size(x, 1),eltype(x)}))
end

@doc doc"""
    log(M::Euclidean, x, y)

Compute the logarithmic map on the [`Euclidean`](@ref) `M` from `x` to `y`,
which in this case is just
````math
\log_x y = y - x.
````
"""
log(::Euclidean, ::Any...)

log!(M::Euclidean, v, x, y) = (v .= y .- x)

log_local_metric_density(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = zero(eltype(x))

@generated _product_of_dimensions(::Euclidean{N}) where {N} = prod(N.parameters)

"""
    manifold_dimension(M::Euclidean)

Return the manifold dimension of the [`Euclidean`](@ref) `M`, i.e.
the product of all array dimensions and the [`real_dimension`](@ref) of the
underlying number system.
"""
function manifold_dimension(M::Euclidean{N,𝔽}) where {N,𝔽}
    return _product_of_dimensions(M) * real_dimension(𝔽)
end

mean(::Euclidean{Tuple{1}}, x::AbstractVector{<:Number}; kwargs...) = mean(x)
function mean(
    ::Euclidean{Tuple{1}},
    x::AbstractVector{<:Number},
    w::AbstractWeights;
    kwargs...,
)
    return mean(x, w)
end
mean(::Euclidean, x::AbstractVector; kwargs...) = mean(x)

function mean!(M::Euclidean, y, x::AbstractVector, w::AbstractVector; kwargs...)
    return mean!(M, y, x, w, GeodesicInterpolation(); kwargs...)
end

function mean_and_var(::Euclidean{Tuple{1}}, x::AbstractVector{<:Number}; kwargs...)
    m, v = mean_and_var(x; kwargs...)
    return m, sum(v)
end
function mean_and_var(
    ::Euclidean{Tuple{1}},
    x::AbstractVector{<:Number},
    w::AbstractWeights;
    corrected = false,
    kwargs...,
)
    m, v = mean_and_var(x, w; corrected = corrected, kwargs...)
    return m, sum(v)
end
function mean_and_var(M::Euclidean, x::AbstractVector, w::AbstractWeights; kwargs...)
    return mean_and_var(M, x, w, GeodesicInterpolation(); kwargs...)
end

median(::Euclidean{Tuple{1}}, x::AbstractVector{<:Number}; kwargs...) = median(x)
function median(
    ::Euclidean{Tuple{1}},
    x::AbstractVector{<:Number},
    w::AbstractWeights;
    kwargs...,
)
    return median(x, w)
end

function median!(::Euclidean{Tuple{1}}, y, x::AbstractVector; kwargs...)
    return copyto!(y, [median(vcat(x...))])
end
function median!(::Euclidean{Tuple{1}}, y, x::AbstractVector, w::AbstractWeights; kwargs...)
    return copyto!(y, [median(vcat(x...), w)])
end

@doc doc"""
    norm(M::Euclidean, x, v)

Compute the norm of a tangent vector `v` at `x` on the [`Euclidean`](@ref)
`M`, i.e. since every tangent space can be identified with `M` itself
in this case, just the (Frobenius) norm of `v`.
"""
norm(::Euclidean, x, v) = norm(v)
norm(::MetricManifold{<:Manifold,EuclideanMetric}, x, v) = norm(v)

"""
    normal_tvector_distribution(M::Euclidean, x, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `x`.
"""
function normal_tvector_distribution(M::Euclidean{Tuple{N}}, x, σ) where {N}
    d = Distributions.MvNormal(zero(x), σ)
    return ProjectedFVectorDistribution(TangentBundleFibers(M), x, d, project_vector!, x)
end

@doc doc"""
    project_point(M::Euclidean, x)

Project an arbitrary point `x` onto the [`Euclidean`](@ref) `M`, which
is of course just the identity map.
"""
project_point(::Euclidean, ::Any...)

project_point!(M::Euclidean, x) = x

"""
    project_tangent(M::Euclidean, x, v)

Project an arbitrary vector `v` into the tangent space of a point `x` on the
[`Euclidean`](@ref) `M`, which is just the identity, since any tangent
space of `M` can be identified with all of `M`.
"""
project_tangent(::Euclidean, ::Any...)

project_tangent!(M::Euclidean, w, x, v) = copyto!(w, v)

"""
    projected_distribution(M::Euclidean, d, [x])

Wrap the standard distribution `d` into a manifold-valued distribution. Generated
points will be of similar type to `x`. By default, the type is not changed.
"""
function projected_distribution(M::Euclidean, d, x)
    return ProjectedPointDistribution(M, d, project_point!, x)
end
function projected_distribution(M::Euclidean, d)
    return ProjectedPointDistribution(M, d, project_point!, rand(d))
end

"""
    representation_size(M::Euclidean)

Return the array dimensions required to represent an element on the
[`Euclidean`](@ref) `M`, i.e. the vector of all array dimensions.
"""
@generated representation_size(::Euclidean{N}) where {N} = size_to_tuple(N)

"""
    sharp(M::Euclidean, x, w)

since cotangent and tangent vectors can directly be identified in the [`Euclidean`](@ref)
case, this yields just the identity for a cotangent vector `w` in the tangent space
of `x` on `M`.
"""
sharp(::Euclidean, ::Any...)

sharp!(M::Euclidean, v::TFVector, x, w::CoTFVector) = copyto!(v, w)

"""
    vector_transport_to(M::Euclidean, x, v, y, ::ParallelTransport)

Parallely transport the vector `v` from the tangent space at `x` to the tangent space at `y`
on the [`Euclidean`](@ref) `M`, which simplifies to the identity.
"""
vector_transport_to(::Euclidean, ::Any, ::Any, ::Any, ::ParallelTransport)

vector_transport_to!(M::Euclidean, vto, x, v, y, ::ParallelTransport) = copyto!(vto, v)

var(::Euclidean, x::AbstractVector; kwargs...) = sum(var(x; kwargs...))
function var(::Euclidean, x::AbstractVector{T}, m::T; kwargs...) where {T}
    return sum(var(x; mean = m, kwargs...))
end

vee(::Euclidean{N,ℝ}, x, v) where {N} = vec(v)

vee!(::Euclidean{N,ℝ}, vⁱ, x, v) where {N} = copyto!(vⁱ, v)

"""
    zero_tangent_vector(M::Euclidean, x)

Return the zero vector in the tangent space of `x` on the [`Euclidean`](@ref)
`M`, which here is just a zero filled array the same size as `x`.
"""
zero_tangent_vector(::Euclidean, ::Any...)

zero_tangent_vector!(M::Euclidean, v, x) = fill!(v, 0)


text:
@inline stridedpointer_for_broadcast(A) = stridedpointer_for_broadcast(ArrayInterface.size(A), stridedpointer(A))
@inline stridedpointer_for_broadcast(s, ptr) = ptr
# function stridedpointer_for_broadcast(s, ptr::VectorizationBase.AbstractStridedPointer)
#      # FIXME: this is unsafe for AbstractStridedPointers
#     throw("Broadcasting not currently supported for arrays where typeof(stridedpointer(A)) === $(typeof(ptr))")
# end
function stridedpointer_for_broadcast_quote(typ, N, S, X)
    q = Expr(:block, Expr(:meta,:inline), :(strd = ptr.strd))
    strd_tup = Expr(:tuple)
    for n ∈ 1:N
        s_type = S[n]
        if s_type <: Static
            if s_type === Static{1}
                push!(strd_tup.args, Expr(:call, lv(:Zero)))
            else
                push!(strd_tup.args, :(strd[$n]))
            end
        else
            Xₙ_type = X[n]
            if Xₙ_type <: Static # FIXME; what to do here? Dynamic dispatch?
                push!(strd_tup.args, :(strd[$n]))
            else
                push!(strd_tup.args, :(Base.ifelse(isone(s[$n]), zero($Xₙ_type), strd[$n])))
            end
        end
    end
    push!(q.args, :(@inbounds $typ(ptr.p, $strd_tup, ptr.offsets)))
    q
end
@generated function stridedpointer_for_broadcast(s::Tuple{Vararg{Any,N}}, ptr::StridedPointer{T,N,C,B,R,X,O}) where {T,N,C,B,R,X,O}
    typ = Expr(:curly, :StridedPointer, T, N, C, B, R)
    stridedpointer_for_broadcast_quote(typ, N, s.parameters, X.parameters)
end
@generated function stridedpointer_for_broadcast(s::Tuple{Vararg{Any,N}}, ptr::VectorizationBase.StridedBitPointer{N,C,B,R,X,O}) where {N,C,B,R,X,O}
    typ = Expr(:curly, lv(:StridedBitPointer), N, C, B, R)
    stridedpointer_for_broadcast_quote(typ, N, s.parameters, X.parameters)
end

struct Product{A,B}
    a::A
    b::B
end
@inline function Base.size(p::Product)
    M = @inbounds size(p.a)[1]
    (M, Base.tail(size(p.b))...)
end
@inline function Base.size(p::Product, i::Integer)
    i == 1 && return @inbounds size(p.a)[1]
    @inbounds size(p.b)[i]
end
@inline Base.length(p::Product) = prod(size(p))
@inline Base.broadcastable(p::Product) = p
@inline numdims(A) = ndims(A) # fallback
@inline numdims(::Type{Product{A,B}}) where {A,B} = numdims(B)
@inline Base.ndims(::Type{Product{A,B}}) where {A,B} = numdims(B)
# This numdims nonsense is a hack to avoid type piracy in defining:
@inline numdims(::Type{B}) where {N, S <: Base.Broadcast.AbstractArrayStyle{N}, B <: Base.Broadcast.Broadcasted{S}} = N

Base.Broadcast._broadcast_getindex_eltype(::Product{A,B}) where {T, A <: AbstractVecOrMat{T}, B <: AbstractVecOrMat{T}} = T
function Base.Broadcast._broadcast_getindex_eltype(p::Product)
    promote_type(
        Base.Broadcast._broadcast_getindex_eltype(p.a),
        Base.Broadcast._broadcast_getindex_eltype(p.b)
    )
end

# recursive_eltype(::Type{A}) where {T, A <: AbstractArray{T}} = T
# recursive_eltype(::Type{NTuple{N,T}}) where {N,T<:Union{Float32,Float64}} = T
# recursive_eltype(::Type{Float32}) = Float32
# recursive_eltype(::Type{Float64}) = Float64
# recursive_eltype(::Type{Tuple{T}}) where {T} = T
# recursive_eltype(::Type{Tuple{T1,T2}}) where {T1,T2} = promote_type(recursive_eltype(T1), recursive_eltype(T2))
# recursive_eltype(::Type{Tuple{T1,T2,T3}}) where {T1,T2,T3} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4}}) where {T1,T2,T3,T4} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4))
# recursive_eltype(::Type{Tuple{T1,T2,T3,T4,T5}}) where {T1,T2,T3,T4,T5} = promote_type(recursive_eltype(T1), recursive_eltype(T2), recursive_eltype(T3), recursive_eltype(T4), recursive_eltype(T5))

# function recursive_eltype(::Type{Broadcasted{S,A,F,ARGS}}) where {S,A,F,ARGS}
#     recursive_eltype(ARGS)
# end

"""
    A *ˡ B

A lazy product of `A` and `B`. While functionally identical to `A * B`, this may avoid the
need for intermediate storage for any computations in `A` or `B`.  Example:

    @avx @. a + B *ˡ (c + d')

which is equivalent to

     a .+ B * (c .+ d')

It should only be used inside an `@avx` block, and to materialize the result it cannot be
the final operation.
"""
@inline *ˡ(a::A, b::B) where {A,B} = Product{A,B}(a, b)
@inline Base.Broadcast.broadcasted(::typeof(*ˡ), a::A, b::B) where {A, B} = Product{A,B}(a, b)
# TODO: Need to make this handle A or B being (1 or 2)-D broadcast objects.
function add_broadcast!(
    ls::LoopSet, mC::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(prod::Type{<:Product}), elementbytes::Int
)
    A, B = prod.parameters
    K = gensym(:K)
    mA = gensym(:Aₘₖ)
    mB = gensym(:Bₖₙ)
    pushprepreamble!(ls, Expr(:(=), mA, Expr(:(.), bcname, QuoteNode(:a))))
    pushprepreamble!(ls, Expr(:(=), mB, Expr(:(.), bcname, QuoteNode(:b))))
    pushprepreamble!(ls, Expr(:(=), K, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), Expr(:ref, Expr(:call, :size, mB), 1))))
    k = gensym(:k)
    add_loop!(ls, Loop(k, 1, K), k)
    m = loopsyms[1];
    if numdims(B) == 1
        bloopsyms = Symbol[k]
        cloopsyms = Symbol[m]
        reductdeps = Symbol[m, k]
        kvec = bloopsyms
    elseif numdims(B) == 2
        n = loopsyms[2];
        bloopsyms = Symbol[k,n]
        cloopsyms = Symbol[m,n]
        reductdeps = Symbol[m, k, n]
        kvec = Symbol[k]
    else
        throw("B must be a vector or matrix.")
    end
    # load A
    # loadA = add_load!(ls, gensym(:A), productref(A, mA, m, k), elementbytes)
    loadA = add_broadcast!(ls, gensym(:A), mA, Symbol[m,k], A, elementbytes)
    # load B
    loadB = add_broadcast!(ls, gensym(:B), mB, bloopsyms, B, elementbytes)
    # set Cₘₙ = 0
    # setC = add_constant!(ls, zero(promote_type(recursive_eltype(A), recursive_eltype(B))), cloopsyms, mC, elementbytes)
    # targetC will be used for reduce_to_add
    mCt = gensym(mC)
    targetC = add_constant!(ls, gensym(:zero), cloopsyms, mCt, elementbytes, :numericconstant)
    push!(ls.preamble_zeros, (identifier(targetC), IntOrFloat))
    setC = add_constant!(ls, gensym(:zero), cloopsyms, mC, elementbytes, :numericconstant)
    push!(ls.preamble_zeros, (identifier(setC), IntOrFloat))
    setC.reduced_children = kvec
    # compute Cₘₙ += Aₘₖ * Bₖₙ
    instrsym = Base.libllvm_version < v"11.0.0" ? :vfmadd231 : :vfmadd
    reductop = Operation(
        ls, mC, elementbytes, instrsym, compute, reductdeps, kvec, Operation[loadA, loadB, setC]
    )
    reductop = pushop!(ls, reductop, mC)
    reductfinal = Operation(
        ls, mCt, elementbytes, :reduce_to_add, compute, cloopsyms, kvec, Operation[reductop, targetC]
    )
    pushop!(ls, reductfinal, mCt)
end

struct LowDimArray{D,T,N,A<:DenseArray{T,N}} <: DenseArray{T,N}
    data::A
end
Base.@propagate_inbounds Base.getindex(A::LowDimArray, i...) = getindex(A.data, i...)
@inline Base.size(A::LowDimArray) = Base.size(A.data)
@inline Base.size(A::LowDimArray, i) = Base.size(A.data, i)
@inline Base.strides(A::LowDimArray) = strides(A.data)
@inline ArrayInterface.parent_type(::Type{LowDimArray{D,T,N,A}}) where {T,D,N,A} = A
@inline ArrayInterface.strides(A::LowDimArray) = ArrayInterface.strides(A.data)
@generated function ArrayInterface.size(A::LowDimArray{D,T,N}) where {D,T,N}
    t = Expr(:tuple)
    for n ∈ 1:N
        if n > length(D) || D[n]
            push!(t.args, Expr(:ref, :s, n))
        else
            push!(t.args, Expr(:call, Expr(:curly, lv(:Static), 1)))
        end
    end
    Expr(:block, Expr(:meta,:inline), :(s = ArrayInterface.size(parent(A))), t)
end
Base.parent(A::LowDimArray) = A.data
Base.unsafe_convert(::Type{Ptr{T}}, A::LowDimArray{D,T}) where {D,T} = pointer(A.data)
ArrayInterface.contiguous_axis(A::LowDimArray) = ArrayInterface.contiguous_axis(A.data)
ArrayInterface.contiguous_batch_size(A::LowDimArray) = ArrayInterface.contiguous_batch_size(A.data)
ArrayInterface.stride_rank(A::LowDimArray) = ArrayInterface.stride_rank(A.data)
ArrayInterface.offsets(A::LowDimArray) = ArrayInterface.offsets(A.data)

@inline function stridedpointer_for_broadcast(A::LowDimArray{D}) where {D}
    _stridedpointer(stridedpointer_for_broadcast(parent(A)), Val{D}())
end

@generated function _stridedpointer(p::StridedPointer{T,N,C,B,R}, ::Val{D}) where {T,N,C,B,R,D}
    lenD = length(D)
    strd = Expr(:tuple)
    offsets = Expr(:tuple)
    Rtup = Expr(:tuple)
    Cnew = -1
    Bnew = -1
    Nnew = 0
    for n ∈ 1:N
        ((n ≤ lenD) && (!D[n])) && continue
        if n == C
            Cnew = n
        end
        if n == B
            Bnew = n
        end
        push!(Rtup.args, R[n])
        push!(offsets.args, Expr(:ref, :offs, n))
        push!(strd.args, Expr(:ref, :strd, n))
        Nnew += 1
    end
    typ = Expr(:curly, :StridedPointer, T, Nnew, Cnew, Bnew, Rtup)
    ptr = Expr(:call, typ, :(pointer(p)), strd, offsets)
    Expr(:block, Expr(:meta,:inline), :(strd = p.strd), :(offs = p.offsets), ptr)
end
# @generated function VectorizationBase.stridedpointer(A::LowDimArray{D,T,N}) where {D,T,N}
#     smul = Expr(:(.), Expr(:(.), :LoopVectorization, QuoteNode(:VectorizationBase)), QuoteNode(:staticmul))
#     multup = Expr(:tuple)
#     for n ∈ D[1]+1:N
#         if length(D) < n
#             push!(multup.args, Expr(:call, :ifelse, :(isone(size(A,$n))), 0, Expr(:ref, :strideA, n)))
#         elseif D[n]
#             push!(multup.args, Expr(:ref, :strideA, n))
#         end
#     end
#     s = Expr(:call, smul, T, multup)
#     f = D[1] ? :PackedStridedPointer : :SparseStridedPointer
#     Expr(:block, Expr(:meta,:inline), Expr(:(=), :strideA, Expr(:call, :strides, Expr(:(.), :A, QuoteNode(:data)))),
#          Expr(:call, Expr(:(.), :VectorizationBase, QuoteNode(f)), Expr(:call, :pointer, :A), s))
# end
function LowDimArray{D}(data::A) where {D,T,N,A <: AbstractArray{T,N}}
    LowDimArray{D,T,N,A}(data)
end
function extract_all_1_array!(ls::LoopSet, bcname::Symbol, N::Int, elementbytes::Int)
    refextract = gensym(bcname)
    ref = Expr(:ref, bcname); foreach(_ -> push!(ref.args, :begin), 1:N)
    pushprepreamble!(ls, Expr(:(=), refextract, ref))
    return add_constant!(ls, refextract, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(LDA::Type{LowDimArray{D,T,N,A}}), elementbytes::Int
) where {D,T,N,A}
    # D,T,N::Int,_ = LDA.parameters
    Dlen = length(D)
    if Dlen == N && !any(D) # array is a scalar, as it is broadcasted on all dimensions
        return extract_all_1_array!(ls, bcname, N, elementbytes)
    end
    fulldims = Symbol[loopsyms[n] for n ∈ 1:N if ((Dlen < n) || D[n]::Bool)]
    ref = ArrayReference(bcname, fulldims)
    add_simple_load!(ls, destname, ref, elementbytes, true, true )::Operation
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{A}, elementbytes::Int
) where {T,N,A<:AbstractArray{T,N}}
    parent = gensym(:parent)
    pushprepreamble!(ls, Expr(:(=), parent, Expr(:call, :parent, bcname)))
    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
    ref = ArrayReference(parent, Symbol[loopsyms[N + 1 - n] for n ∈ 1:N])
    add_simple_load!( ls, destname, ref, elementbytes, true, true )::Operation
end
function add_broadcast_adjoint_array!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{<:AbstractVector}, elementbytes::Int
)
    # isone(length(loopsyms)) && return extract_all_1_array!(ls, bcname, N, elementbytes)
    parent = gensym(:parent)
    pushprepreamble!(ls, Expr(:(=), parent, Expr(:call, :parent, bcname)))

    ref = ArrayReference(parent, Symbol[loopsyms[2]])
    add_simple_load!( ls, destname, ref, elementbytes, true, true )::Operation
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Adjoint{T,A}}, elementbytes::Int
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{Transpose{T,A}}, elementbytes::Int
) where {T, A <: AbstractArray{T}}
    add_broadcast_adjoint_array!( ls, destname, bcname, loopsyms, A, elementbytes )
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    ::Type{<:AbstractArray{T,N}}, elementbytes::Int
) where {T,N}
    add_simple_load!(ls, destname, ArrayReference(bcname, @view(loopsyms[1:N])), elementbytes, true, true)
end
function add_broadcast!(
    ls::LoopSet, ::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{T}, elementbytes::Int
) where {T<:Number}
    add_constant!(ls, bcname, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, ::Symbol, bcname::Symbol, loopsyms::Vector{Symbol}, ::Type{Base.RefValue{T}}, elementbytes::Int
) where {T}
    refextract = gensym(bcname)
    pushprepreamble!(ls, Expr(:(=), refextract, Expr(:ref, bcname)))
    add_constant!(ls, refextract, elementbytes) # or replace elementbytes with sizeof(T) ? u
end
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(_::Type{SubArray{T,N,A,S,B}}), elementbytes::Int
) where {T,N,N2,A<:AbstractArray{T,N2},B,N3,S <: Tuple{Int,Vararg{Any,N3}}}
    inds = Vector{Symbol}(undef, N+1)
    inds[1] = Symbol("##DISCONTIGUOUSSUBARRAY##")
    inds[2:end] .= @view(loopsyms[1:N])
    add_simple_load!(ls, destname, ArrayReference(bcname, inds), elementbytes, true, true)
end
const BroadcastedArray{S<:Broadcast.AbstractArrayStyle,F,A} = Broadcasted{S,Nothing,F,A}
function add_broadcast!(
    ls::LoopSet, destname::Symbol, bcname::Symbol, loopsyms::Vector{Symbol},
    @nospecialize(B::Type{<:BroadcastedArray}),
    elementbytes::Int
)
    S,_,F,A = B.parameters
    instr = get(FUNCTIONSYMBOLS, F) do
        f = gensym(:func)
        pushprepreamble!(ls, Expr(:(=), f, Expr(:(.), bcname, QuoteNode(:f))))
        Instruction(bcname, f)
    end
    args = A.parameters
    Nargs = length(args)
    bcargs = Expr(:(.), bcname, QuoteNode(:args))
    # this is the var name in the loop
    parents = Operation[]
    deps = Symbol[]
    # reduceddeps = Symbol[]
    for (i,arg) ∈ enumerate(args)
        argname = gensym(:arg)
        pushprepreamble!(ls, Expr(:(=), argname, Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), Expr(:ref, bcargs, i))))
        # dynamic dispatch
        parent = add_broadcast!(ls, gensym(:temp), argname, loopsyms, arg, elementbytes)::Operation
        push!(parents, parent)
        mergesetdiffv!(deps, loopdependencies(parent), reduceddependencies(parent))
    end
    op = Operation(
        length(operations(ls)), destname, elementbytes, instr, compute, deps, NOPARENTS, parents
    )
    pushop!(ls, op, destname)
end

# size of dest determines loops
# function vmaterialize!(
@generated function vmaterialize!(
    dest::AbstractArray{T,N}, bc::BC, ::Val{Mod}
) where {T <: NativeTypes, N, BC <: Union{Broadcasted,Product}, Mod}
    # we have an N dimensional loop.
    # need to construct the LoopSet
    # @show typeof(dest)
    loopsyms = [gensym(:n) for n ∈ 1:N]
    ls = LoopSet(Mod)
    ls.isbroadcast[] = true
    sizes = Expr(:tuple)
    for (n,itersym) ∈ enumerate(loopsyms)
        Nsym = gensym(:N)
        add_loop!(ls, Loop(itersym, 1, Nsym), itersym)
        push!(sizes.args, Nsym)
    end
    pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest)))
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    add_simple_store!(ls, :dest, ArrayReference(:dest, loopsyms), elementbytes)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    # return ls
    q = lower(ls, 0)
    push!(q.args, :dest)
    # @show q
    # q
    q = Expr(
        :block,
        ls.prepreamble,
        # Expr(:if, check_args_call(ls), Expr(:block, :(println("Primary code path!")), q), Expr(:block, :(println("Back up code path!")), :(Base.Broadcast.materialize!(dest, bc))))
        Expr(:if, check_args_call(ls), q, :(Base.Broadcast.materialize!(dest, bc)))
    )
    # isone(N) && pushfirst!(q.args, Expr(:meta,:inline))
    q
     # ls
end
@generated function vmaterialize!(
    dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::BC, ::Val{Mod}
) where {T <: NativeTypes, N, A <: AbstractArray{T,N}, BC <: Union{Broadcasted,Product}, Mod}
    # we have an N dimensional loop.
    # need to construct the LoopSet
    loopsyms = [gensym(:n) for n ∈ 1:N]
    ls = LoopSet(Mod)
    ls.isbroadcast[] = true
    pushprepreamble!(ls, Expr(:(=), :dest, Expr(:call, :parent, :dest′)))
    sizes = Expr(:tuple)
    for (n,itersym) ∈ enumerate(loopsyms)
        Nsym = gensym(:N)
        add_loop!(ls, Loop(itersym, 1, Nsym), itersym)
        push!(sizes.args, Nsym)
    end
    pushpreamble!(ls, Expr(:(=), sizes, Expr(:call, :size, :dest′)))
    elementbytes = sizeof(T)
    add_broadcast!(ls, :dest, :bc, loopsyms, BC, elementbytes)
    add_simple_store!(ls, :dest, ArrayReference(:dest, reverse(loopsyms)), elementbytes)
    resize!(ls.loop_order, num_loops(ls)) # num_loops may be greater than N, eg Product
    q = lower(ls, 0)
    push!(q.args, :dest′)
    q = Expr(
        :block,
        ls.prepreamble,
        Expr(:if, check_args_call(ls), q, :(Base.Broadcast.materialize!(dest′, bc)))
    )
    # isone(N) && pushfirst!(q.args, Expr(:meta,:inline))
    q
    # ls
end
# these are marked `@inline` so the `@avx` itself can choose whether or not to inline.
@inline function vmaterialize!(
    dest::AbstractArray{T,N}, bc::Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(identity),Tuple{T2}}, ::Val{Mod}
) where {T <: NativeTypes, N, T2 <: Number, Mod}
    arg = T(first(bc.args))
    @avx for i ∈ eachindex(dest)
        dest[i] = arg
    end
    dest
end
@inline function vmaterialize!(
    dest′::Union{Adjoint{T,A},Transpose{T,A}}, bc::Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(identity),Tuple{T2}}, ::Val{Mod}
) where {T <: NativeTypes, N, A <: AbstractArray{T,N}, T2 <: NativeTypes, Mod}
    arg = T(first(bc.args))
    dest = parent(dest′)
    @avx for i ∈ eachindex(dest)
        dest[i] = arg
    end
    dest′
end

@inline function vmaterialize(bc::Broadcasted, ::Val{Mod}) where {Mod}
    ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
    vmaterialize!(similar(bc, ElType), bc, Val{Mod}())
end

vmaterialize!(dest, bc, ::Val{mod}) where {mod} = Base.Broadcast.materialize!(dest, bc)


text:
module Functions
export sigmoid, relu, softmax, cross_entropy_error

# TODO:関数の内部にdotを付けるのではなく，関数の呼び出し元にdotを付けるように変更
function sigmoid(x)
    return 1 ./ (1 .+ exp.(-x))
end

function relu(x)
    return max(0, x) # NOT maximum()
end

function softmax(x)
    y = []
    for row_i = 1:size(x, 1)
        c = maximum(x[row_i, :])
        exp_a = exp.(x[row_i, :] .- c) # オーバフロー対策
        push!(y, exp_a ./ sum(exp_a))
    end
    return hcat(y...)' # cf. https://docs.julialang.org/en/v1/manual/faq/?highlight=splat#What-does-the-...-operator-do?
end

function cross_entropy_error(y, t)
    return -sum(t .* log.(y .+ 1e-7)) # 0除算を回避するために微小な値を加算
end
end