text:
@testset "Root counts" begin

    @testset "Benchmark systems" begin
        kat10 = solve(equations(katsura(10)); show_progress = false)
        @test nfinite(kat10) == 1024
        @test nsingular(kat10) == 0
        @test nreal(kat10) == 216

        cyclic5 = solve(equations(cyclic(5)); show_progress = false)
        @test nfinite(cyclic5) == 70
        @test nreal(cyclic5) == 10
        @test nsingular(cyclic5) == 0

        cyclic7 = solve(equations(cyclic(7)); show_progress = false)
        @test nfinite(cyclic7) == 924
        @test nreal(cyclic7) == 56
        @test nsingular(cyclic7) == 0
        @test ntracked(cyclic7) == 5040
    end

    @testset "Multi-homogeneous" begin
        ## initialize the variables
        @polyvar z[1:6, 1:3]
        F = let
            p = [1, 1, 0]
            α = [
                0.9606655201575953,
                -0.12634585123792536,
                -0.821994270849398,
                -0.4612529209681106,
                0.1914597581898848,
            ]
            a = [
                -0.6508165773155992,
                -0.6252070694344862,
                1.2753203087262615,
                -0.30371050526688304,
                0.45398822665915634,
                -0.4042256983795828,
                -1.1485430048515646,
                -0.22145594990414993,
                0.9590486265536351,
            ]
            ## define the system of polynomials
            f = [z[i, :] ⋅ z[i, :] for i = 2:5]
            g = [z[i, :] ⋅ z[i+1, :] for i = 1:5]
            h = sum(a[i] .* (z[i, :] × z[i+1, :]) for i = 1:3) +
                sum(a[i+4] .* z[i, :] for i = 2:5)
            F′ = [f .- 1; g .- cos.(α); h .- p]
            ## assign values to z₁ and z₆
            [subs(f, z[1, :] => [1, 0, 0], z[6, :] => [1, 0, 0]) for f in F′]
        end
        group1 = [[z[2, :]; z[4, :]], [z[3, :]; z[5, :]]]
        @test nnonsingular(solve(
            F;
            system = FPSystem,
            seed = 141232,
            show_progress = false,
        )) == 16
        @test nnonsingular(solve(
            F;
            variable_groups = group1,
            system = FPSystem,
            seed = 12312,
            show_progress = false,
        )) == 16
    end

    @testset "Bio-chemical reaction networks" begin
        @testset "bacillus_subtilis" begin
            f = equations(PolynomialTestSystems.bacillus_subtilis())
            res = solve(f; show_progress = false, seed = 173221)
            @test nsolutions(res) == 44
            @test nsingular(res) == 0
            res = solve(f; seed = 231231, start_system = :polyhedral, show_progress = false)
            @test nsolutions(res) == 44
        end

        # These test got prepared by Torkel Loman to test the BioDiffEq integration
        @polyvar x[1:3] p[1:10]
        results1_direct = fill(0, 12)
        results1_template = fill(0, 12)
        results2_direct = fill(0, 12)
        results2_template = fill(0, 12)
        results3_direct = fill(0, 12)
        results3_template = fill(0, 12)
        results4_direct = fill(0, 12)
        results4_template = fill(0, 12)

        for i = 1:100
            pol = [
                -x[1] * x[3] * p[3] - x[1] * p[2] + p[1],
                x[1] * x[2] * x[3] * p[8] * p[9] + x[1] * x[3] * p[7] * p[8] * p[9] -
                x[2] * x[3] * p[5] * p[6] - x[2] * p[5] * p[6] * p[10] -
                x[2] * x[3] * p[4] - x[2] * p[4] * p[10],
                -x[1] * x[2] * x[3] * p[8] * p[9] - x[1] * x[3] * p[7] * p[8] * p[9] +
                x[2] * x[3] * p[5] * p[6] + x[2] * p[5] * p[6] * p[10] +
                x[2] * x[3] * p[4] + x[2] * p[4] * p[10],
                x[2] + x[3] - 1.0,
            ]
            p_vals = [0.04, 0.04, 1.0, 1.0, 10.0, 0.0, 0.04, 35.0, 0.1, 0.04]
            pol_pars = DynamicPolynomials.subs.(pol, Ref(p => p_vals))
            sol = solutions(HomotopyContinuation.solve(pol_pars, show_progress = false))
            results1_direct[length(sol)+1] += 1

            p_template = randn(ComplexF64, 10)
            f_template = DynamicPolynomials.subs.(pol, Ref(p => p_template))
            result_template = solutions(HomotopyContinuation.solve(
                f_template,
                show_progress = false,
            ))
            (length(result_template) > 0) && (sol_again = solutions(HomotopyContinuation.solve(
                pol,
                result_template,
                parameters = p,
                p₁ = p_template,
                p₀ = ComplexF64.(p_vals),
                show_progress = false,
            )))
            (length(result_template) > 0) ? (results1_template[length(sol_again)+1] += 1) :
            (results1_template[1] += 1)

            pol2 = [
                -x[1]^5 * x[2] * p[5] + x[1]^4 * x[3] * p[6] * p[8] -
                x[1] * x[2] * p[3]^4 * p[5] + x[3] * p[3]^4 * p[6] * p[8] - x[1]^5 * p[7] +
                x[1]^4 * x[3] * p[4] - x[1] * p[3]^4 * p[7] +
                x[3] * p[3]^4 * p[4] + x[1]^4 * p[1] + x[1]^4 * p[2] + p[1] * p[3]^4,
                -x[1]^5 * x[2] * p[5] - x[1] * x[2] * p[3]^4 * p[5] - x[1]^4 * x[2] * p[7] +
                x[1]^4 * x[3] * p[4] - x[2] * p[3]^4 * p[7] +
                x[3] * p[3]^4 * p[4] + x[1]^4 * p[1] + x[1]^4 * p[2] + p[1] * p[3]^4,
                x[1] * x[2] * p[5] - x[3] * p[6] * p[8] - x[3] * p[4] - x[3] * p[7],
            ]
            p2_vals = [0.005, 0.1, 2.8, 10, 100, 0.1, 0.01, 0.0]
            pol2_pars = DynamicPolynomials.subs.(pol2, Ref(p => p2_vals))
            sol2 = solutions(HomotopyContinuation.solve(pol2_pars, show_progress = false))
            results2_direct[length(sol2)+1] += 1

            p2_template = randn(ComplexF64, 8)
            f2_template = DynamicPolynomials.subs.(pol2, Ref(p => p2_template))
            solve_templ2 = HomotopyContinuation.solve(f2_template, show_progress = false)
            result2_template = solutions(solve_templ2)
            sol2_again = solutions(HomotopyContinuation.solve(
                pol2,
                result2_template,
                parameters = p[1:8],
                p₁ = p2_template,
                p₀ = ComplexF64.(p2_vals),
                precision_strategy = :adaptive_always,
                max_steps = 10_000,
                show_progress = false,
            ))
            results2_template[length(sol2_again)+1] += 1

            pol3 = pol2
            p3_vals = [0.005, 0.1, 2.8, 10, 100, 0.1, 0.01, 1.0] #Last parameter is 1 and not 0. Here there should be 3 real solutions instead of 1.
            pol3_pars = DynamicPolynomials.subs.(pol3, Ref(p => p3_vals))
            sol3 = solutions(HomotopyContinuation.solve(pol3_pars, show_progress = false))
            results3_direct[length(sol3)+1] += 1

            p3_template = randn(ComplexF64, 8)
            f3_template = DynamicPolynomials.subs.(pol3, Ref(p => p3_template))
            result3_template = solutions(HomotopyContinuation.solve(
                f3_template,
                show_progress = false,
            ))
            sol3_again = solutions(HomotopyContinuation.solve(
                pol3,
                result3_template,
                parameters = p[1:8],
                p₁ = p3_template,
                p₀ = ComplexF64.(p3_vals),
                precision_strategy = :adaptive_always,
                max_steps = 10_000,
                show_progress = false,
            ))
            results3_template[length(sol3_again)+1] += 1

            pol4 = [-x[1]^3 * p[3] - x[1] * p[2]^2 * p[3] + x[1]^2 * p[1]]
            p4_vals = [1.0, 0.2, 1.0]
            pol4_pars = DynamicPolynomials.subs.(pol4, Ref(p => p4_vals))
            sol4 = solutions(HomotopyContinuation.solve(pol4_pars, show_progress = false))
            results4_direct[length(sol4)+1] += 1

            p4_template = randn(ComplexF64, 3)
            f4_template = DynamicPolynomials.subs.(pol4, Ref(p => p4_template))
            result4_template = solutions(HomotopyContinuation.solve(
                f4_template,
                show_progress = false,
            ))
            sol4_again = solutions(HomotopyContinuation.solve(
                pol4,
                result4_template,
                parameters = p[1:3],
                p₁ = p4_template,
                p₀ = ComplexF64.(p4_vals),
                precision_strategy = :adaptive_always,
                max_steps = 10_000,
                show_progress = false,
            ))
            results4_template[length(sol4_again)+1] += 1
        end

        @test 97 ≤ results1_template[4] ≤ 100
        @test 97 ≤ results1_direct[4] ≤ 100
        @test 97 ≤ results2_direct[7] ≤ 100
        @test 97 ≤ results2_template[7] ≤ 100
        @test 97 ≤ results3_direct[7] ≤ 100
        @test 97 ≤ results3_template[7] ≤ 100
        @test 97 ≤ results4_direct[4] ≤ 100
        @test 97 ≤ results4_template[4] ≤ 100
    end
end


text:
data = [1 2; 1 2; 5 7; 2 3]
bar(data)

ext: jl


text:
using Unitful

import Unitful:
    nm, μm, mm, cm, m, km,
    mg, g, kg,
    ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
    eV,
    μJ, mJ, J,
	μW, mW, W

	"""
	GaussianLaser

	Representation of a Gaussian laser defined by a laser, the location of the waist (z0) and the waist radius (w0)

	# Fields
	- `laser::Laser`       : A Laser
	- `w0::Unitful.Length` : radius of the waist
	- `z0::Unitful.Length` : location of the waist (computed from w0)
	- `I0::typeof(1.0mW/cm^2)` : Intensity of the beam at (0,0) (computed from w0)
	- `θ0::Real ` : divergence of the beam  (computed from w0)

	"""
	struct GaussianLaser
		laser::Laser
		w0::Unitful.Length
		z0::Unitful.Length
		I0::typeof(1.0mW/cm^2)
		γ0::typeof(1.0Hz/cm^2)
		θ0::Real

		function GaussianLaser(laser, w0)
			z0 = π * w0^2 / laser.λ
			I0 = 2 * laser.P / (π * w0^2)
			nγ = uconvert(Hz, laser.P / photon_energy(laser.λ))
			γ0 = 2 * nγ / (π * w0^2)
			θ0 = w0/z0
			new(laser, w0, z0, I0, γ0, θ0)
		end
	end

	"""
	W(gl::GaussianLaser)

returns the beam width: ``W(z) = W_0 \\sqrt{1 + (z/z_0)^2}``

# Fields
	- `gl::GaussianLaser`       : A gaussian Laser

"""
function W(gl::GaussianLaser)
	Wz(z::Unitful.Length) = gl.w0 * sqrt(1.0 + (z/gl.z0)^2)
	return Wz
end

"""
R(gl::GaussianLaser)

returns the beam radius of curvature: ``R(z) = z ( 1 + (z_0/z)^2) ``

# Fields
- `gl::GaussianLaser`       : A gaussian Laser

"""
function R(gl::GaussianLaser)
	Rz(z::Unitful.Length) = z * (1.0 + (gl.z0/z)^2)
	return Rz
end

"""
I(gl::GaussianLaser)

returns the beam Intensity: ``I(\\rho, z) = I_0 ( W_0 / W(z))^2 \\exp{-2 \\rho^2/W^2(z)}  ``

# Fields
- `gl::GaussianLaser`       : A gaussian Laser

"""
function I(gl::GaussianLaser)
	Wz = W(gl)
	Irz(ρ::Unitful.Length, z::Unitful.Length) = gl.I0 * (gl.w0 / Wz(z))^2 * exp(-2.0 * ρ^2/Wz(z)^2)
	return Irz
end

"""
	spot_size(gl::GaussianLaser)

returns 2*w0 that is the diameter of the beam waist, also called spot size

# Fields

- `gl::GaussianLaser` : A gaussian laser

"""
spot_size(gl::GaussianLaser)  = 2 * gl.w0

"""
	angular_divergence(gl::GaussianLaser)

returns 2θ that is the beam angular divergence

# Fields

- `gl::GaussianLaser` : A gaussian laser

"""
angular_divergence(gl::GaussianLaser)  = 2 * gl.θ0

"""
	depth_of_focus(gl::GaussianLaser)

returns 2z0

# Fields

- `gl::GaussianLaser` : A gaussian laser

"""
depth_of_focus(gl::GaussianLaser)  = 2 * gl.z0

"""
	propagate_paralel_beam(laser::Laser, obj::Objective)

returns a gaussian laser with parameter w0 resulting from focusing the paralell laser beam
filling the entrance pupil of the objective.

# Fields

- `laser::Laser` : A  laser filling the entrance pupil of the objective (assuming parallell rays)
- `obj::Objective` : An objective, used to focuse the laser in the sample.

"""
function propagate_paralell_beam(laser::Laser, obj::Objective)
	w0 = laser.λ / (π * obj.NA)
	GaussianLaser(laser, w0)
end

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
module MyExample

export f

f(x,y) = 2x + 3y
end


text:
import Base: @deprecate

## deprecated in v0.13

## deprecated in v0.12

@deprecate(Cart_get(comm, maxdims), Cart_get(comm), false)

@deprecate(Cart_coords!(comm, rank, maxdims::Integer, coords),
           Cart_coords!(comm, rank, coords), false)

@deprecate(Cart_coords(comm, maxdims::Integer), Cart_coords(comm), false)


text:
import Random
import ForwardDiff
import DiffResults
using LinearAlgebra
using DataFrames
using CSV
using Statistics
using Printf
include("conjutilities.jl")
# include("grad_hess.jl")
Random.seed!(1)
num_samples = 10
cone_ds = [20, 40, 60]
offsets = 10.0 .^ collect(-5:-1)

##
abstract type Cone end
struct SumLog <: Cone
    d::Int
end
struct HypoPower <: Cone
    α::Vector{Float64}
end
HypoPower(d::Int) = (α = rand(d); α /= sum(α); HypoPower(α))
struct HypoGeom <: Cone
    d::Int
end
struct RadialPower <: Cone
    α::Vector{Float64}
end
RadialPower(d::Int) = (α = rand(d); α /= sum(α); RadialPower(α))
struct RadialGeom <: Cone
    d::Int
end
struct InfinityNorm <: Cone
    d::Int
end

##
function barrier(::SumLog)
    function bar(uvw)
        (u, v, w) = (uvw[1], uvw[2], uvw[3:end])
        return -log(v * sum(log(wi / v) for wi in w) - u) - sum(log, w) - log(v)
    end
    return bar
end
function barrier(cone::HypoPower)
    α = cone.α
    function bar(uw)
        (u, w) = (uw[1], uw[2:end])
        return -log(exp(sum(αi * log(wi) for (αi, wi) in zip(α, w))) - u) -
            sum(log, w)
    end
    return bar
end
function barrier(cone::HypoGeom)
    function bar(uw)
        (u, w) = (uw[1], uw[2:end])
        return -log(exp(sum(log(wi) / cone.d for wi in w)) - u) - sum(log, w)
    end
    return bar
end
function barrier(cone::RadialPower)
    α = cone.α
    function bar(uw)
        (u, w) = (uw[1], uw[2:end])
        return -log(exp(sum(2αi * log(wi) for (αi, wi) in zip(α, w))) - u^2) -
            sum((1 - αi) * log(wi) for (αi, wi) in zip(α, w))
    end
    return bar
end
function barrier(cone::RadialGeom)
    function bar(uw)
        (u, w) = (uw[1], uw[2:end])
        d = cone.d
        return -log(exp(sum(2 / d * log(wi) for wi in w)) - u^2) - (d - 1) / d *
            sum(log(wi) for wi in w)
    end
    return bar
end
function barrier(cone::InfinityNorm)
    function bar(uw)
        (u, w) = (uw[1], uw[2:end])
        return sum(-log(u^2 - wi^2) for wi in w) + (cone.d - 1) * log(u)
    end
    return bar
end

init_guess(cone::SumLog) = vcat(-1, 1, ones(cone.d))
init_guess(cone::HypoPower) = vcat(-1, ones(length(cone.α)))
init_guess(cone::HypoGeom) = vcat(-1, ones(cone.d))
init_guess(cone::RadialPower) = vcat(0, sqrt.(1 .+ cone.α))
init_guess(cone::RadialGeom) = vcat(0, fill(sqrt.(1 + inv(cone.d)), cone.d))
init_guess(cone::InfinityNorm) = vcat(1, zeros(cone.d))

function dual_point(cone::SumLog, offset::Float64)
    d = cone.d
    r = rand(d)
    p = -rand()
    phi = (p * sum(log(-ri / p) for ri in r) + p * d)
    q = phi * (1 + sign(phi) * offset)
    return vcat(p, q, r)
end
function dual_point(cone::HypoPower, offset::Float64)
    α = cone.α
    r = rand(length(α))
    phi = -exp(sum(αi * log(ri / αi) for (ri, αi) in zip(r, α)))
    p = phi * (1 - offset)
    return vcat(p, r)
end
dual_point(cone::HypoGeom, offset::Float64) = dual_point(HypoPower(cone.d), offset)
function dual_point(cone::RadialPower, offset::Float64)
    α = cone.α
    r = rand(length(α))
    # only look at the p > 0 case, but that's OK because residuals aren't
    # going to behave differently in the other direction
    phi = exp(sum(αi * log(ri / αi) for (ri, αi) in zip(r, α)))
    p = phi * (1 - offset)
    return vcat(p, r)
end
dual_point(cone::RadialGeom, offset::Float64) = dual_point(RadialPower(cone.d), offset)
function dual_point(cone::InfinityNorm, offset::Float64)
    r = randn(cone.d)
    p = sum(abs, r) * (1 + offset)
    return vcat(p, r)
end

function grad_hess(cone::Cone, s)
    bar = barrier(cone)
    # s = BigFloat.(s)
    result = DiffResults.HessianResult(s)
    return ForwardDiff.hessian!(result, bar, s);
end

function dual_grad(cone::SumLog, pqr)
    d = cone.d
    (p, q, r) = (pqr[1], pqr[2], pqr[3:end])
    dual_ϕ = sum(log(r_i / -p) for r_i in r)
    β = 1 + d - q / p + dual_ϕ
    bomega = d * wrightomega(β / d - log(d))
    @assert bomega + d * log(bomega) ≈ β

    cgp = (-d - 2 + q / p + 2 * bomega) / (p * (1 - bomega))
    cgq = -1 / (p * (1 - bomega))
    cgr = bomega ./ r / (1 - bomega)
    return (vcat(cgp, cgq, cgr), 0)
end
function dual_grad(cone::HypoPower, pr)
    (p, r) = (pr[1], pr[2:end])
    α = cone.α
    d = length(α)
    sumlog = sum(α_i * log(r_i) for (α_i, r_i) in zip(α, r))

    h(y) = sum(αi * log(y - p * αi) for αi in α) - sumlog
    hp(y) = sum(αi / (y - p * αi) for αi in α)
    lower = 0.0
    upper = exp(sumlog) + p / d
    # (new_bound, iter) = rootnewton(h, hp, lower = lower, upper = upper)
    (new_bound, iter) = rootnewton(h, hp, init = 0.0, increasing = true)

    dual_g_ϕ = inv(new_bound)
    cgp = -inv(p) - dual_g_ϕ
    cgr = (p * dual_g_ϕ * α .- 1) ./ r
    return (vcat(cgp, cgr), iter)
end
function dual_grad(cone::HypoGeom, pr)
    d = cone.d
    (p, r) = (pr[1], pr[2:end])
    dual_ϕ = exp(sum(log, r) / d)
    dual_ζ = dual_ϕ + p / d
    cgp = -inv(p) - inv(dual_ζ)
    cgr = -dual_ϕ / dual_ζ ./ r
    return (vcat(cgp, cgr), 0)
end
function dual_grad(cone::RadialPower, pr)
    (p, r) = (pr[1], pr[2:end])
    α = cone.α
    d = length(α)
    if iszero(p)
        cgp = 0
        cgr = -(α .+ 1) ./ r
    else
        log_phi_r = 2 * sum(αi * log(ri) for (αi, ri) in zip(α, r))
        phi_r = exp(log_phi_r)
        phi_αr = exp(sum(αi * log(ri / αi) for (αi, ri) in zip(α, r)))
        inner_bound = -1 / p - (1 + sign(p) * 1 / p * sqrt(phi_r * (d^2 / p^2 *
            phi_r + d^2 - 1))) / (p / d - d * phi_r / p)
        gamma = abs(p) / phi_αr
        outer_bound = (1 + d) * gamma / (1 - gamma) / p
        h(y) = 2 * sum(αi * log(2 * αi * y^2 + (1 + αi) * 2 * y / p) for αi in α) -
            log_phi_r - log(2 * y / p + y^2) - 2 * log(2 * y / p)
        hp(y) = 2 * sum(αi^2 / (αi * y + (1 + αi) / p) for αi in α) -
            2 * (y + 1 / p) / y / (y + 2 / p)
        (cgp, iter) = rootnewton(h, hp, lower = inner_bound, upper = outer_bound,
            init = outer_bound, increasing = false)
        cgr = -(α * (1 + p * cgp) .+ 1) ./ r
    end
    return (vcat(cgp, cgr), iter)
end
function dual_grad(cone::RadialGeom, pr)
    d = cone.d
    (p, r) = (pr[1], pr[2:end])
    if iszero(p)
        ctp = 0
        cgr = -(inv(d) + 1) ./ r
    else
        phi_r = exp(2 * sum(log(ri) / d for ri in r))
        cgp = -1 / p - (1 + sign(p) * 1 / p * sqrt(phi_r * ((d / p)^2 * phi_r +
            d^2 - 1))) / (p / d - d * phi_r / p)
        cgr = -((1 + p * cgp) / d + 1) ./ r
    end
    return (vcat(cgp, cgr), 0)
end
function dual_grad(cone::InfinityNorm, pr)
    d = cone.d
    (p, r) = (pr[1], pr[2:end])
    dual_zeta = p - sum(abs, r)
    h(y) = p * y + sum(sqrt(1 + abs2(ri * y)) for ri in r) + 1
    hp(y) = p + y * sum(abs2(ri) / sqrt(1 + abs2(ri * y)) for ri in r)
    # lower = -(d + 1) / dual_zeta
    upper = min(-inv(dual_zeta), -(d + 1) / p)
    (cgp, iter) = rootnewton(h, hp, init = upper)

    cgr = copy(r)
    for i in eachindex(r)
        if abs(r[i]) .< eps()
            cgr[i] = 0
        else
            cgr[i] = (-1 + sqrt(1 + abs2(r[i] * cgp))) / r[i]
        end
    end
    return (vcat(cgp, cgr), iter)
end

function naive_dual_grad(cone, z)
    curr = init_guess(cone)
    derivs = grad_hess(cone, curr)
    quad_bound = 0.35
    (g, H) = (DiffResults.gradient(derivs), DiffResults.hessian(derivs))
    H_fact = cholesky(Symmetric(H))
    r = z + g
    Hir = H_fact \ r
    n = n_prev = sqrt(dot(r, Hir))
    #
    # Hiz = H_fact \ z
    # Hir = Hiz - curr
    # n = n_prev = sqrt(length(z) - 2 * dot(curr, z) + dot(z, Hiz))
    max_iter = 400
    iter = 1
    resid_path = []

    while n > 1000eps()
        α = (n > quad_bound ? inv(1 + n) : 1)
        curr -= Hir * α
        push!(resid_path, dot(z, -curr) + length(z))
        derivs = grad_hess(cone, curr)
        (g, H) = (DiffResults.gradient(derivs), DiffResults.hessian(derivs))
        H_fact = cholesky(Symmetric(H), check = false)
        if !issuccess(H_fact)
            H_fact = bunchkaufman(Symmetric(H))
        end
        r = z + g
        Hir = H_fact \ r
        n2 = dot(r, Hir)
        # Hiz = H_fact \ z
        # Hir = Hiz - curr
        # n2 = length(z) - 2 * dot(curr, z) + dot(z, Hiz)
        n2 < -sqrt(eps()) && break # couple crazy hypogeom cases
        n = sqrt(abs(n2))
        iter += 1
        # slow progress (numerical)
        if (n_prev < quad_bound) && (n > 1000(n_prev / (1 - n_prev))^2)
            break
        end
        n_prev = n
        (iter == max_iter) && break
    end
    return (-curr, iter, resid_path)
end

##
cone_ts = [
    SumLog
    HypoPower
    HypoGeom
    RadialPower
    RadialGeom
    InfinityNorm
    ]

function printx(x::Float64)
    if abs(x) < 10
        return @sprintf("%.1f", x)
    else
        return @sprintf("%.0f.", x)
    end
end

function get_resids()
    reio = open("tex/resid.tex", "w")
    itio = open("tex/iters.tex", "w")
    sep = " & "
    for d in cone_ds
        results = DataFrame(
            cone_d = Int[],
            cone_type = Symbol[],
            offset = Float64[],
            trial = Int[],
            residdirect = Float64[],
            residnaive = Float64[],
            newtoniters = Int[],
            directiters = Int[],
            )

        for C in cone_ts, o in offsets, i in 1:num_samples
            Random.seed!(i)
            cone = C(d)
            z = dual_point(cone, o)
            (g, iter) = dual_grad(cone, z)
            resid = abs(dot(z, g) + length(z)) # if -g slightly out of primal cone then can be numerically negative
            # @show resid
            (naive_g, naive_iter, resid_path) = naive_dual_grad(cone, z)
            naive_resid = abs(dot(z, naive_g) + length(z))
            push!(results, (d, nameof(C), o, i, resid,
                naive_resid, naive_iter, iter))
            if i == 1 && d == 60
                @show nameof(C), resid
                resid_df = DataFrame(resid = resid_path)
                resids_name = lowercase(string(nameof(C))) * string(Int(-log10(o)))
                CSV.write(joinpath("csvs", resids_name * ".csv"), resid_df)
            end
        end

        agg = combine(groupby(results, [:cone_d, :cone_type, :offset]),
            :residdirect => mean => :residdirectmean,
            :residnaive => mean => :residnaivemean,
            :newtoniters => mean => :newtonitersmean,
            :directiters => mean => :directitersmean,
            )

        for subdf in groupby(agg, [:cone_d, :offset])
            for r in eachrow(subdf)
                print(reio, sep * printx(log10(abs(r.residnaivemean))) * sep *
                    printx(log10(abs(r.residdirectmean))))
                print(itio, sep * printx(r.newtonitersmean) *
                    sep * printx(r.directitersmean))
            end
            print(reio, "\\\\ \n")
            print(itio, "\\\\ \n")
        end

        CSV.write("csvs/cgs_$(d).csv", agg)
    end
    close(reio)
    close(itio)
    return
end

get_resids()

;

max_stars_repo_name: lkapelevich/conjgradexperiments

max_issues_repo_name: lkapelevich/conjgradexperiments

max_forks_repo_name: lkapelevich/conjgradexperiments


text:
module ModelDerivatives

using ForwardDiff
using Morbit
using Test

using Combinatorics

const scalar_test_funcs = [
	x -> sum(x),
	x -> sum( x.^2 ),
	x -> sum( x.^(2 : (length(x) + 1)) ),
	x -> exp( sum( x.^2 ) / 100 ),
	x -> ( 1 / (10 + abs(sum(x))) ) * sin( sum(x) )
]

const configs = [
	ExactConfig(;gradients = :autodiff),
	ExactConfig(;gradients = :fdm),
	LagrangeConfig(; degree = 1 ),
	LagrangeConfig(; degree = 2 ),
	RbfConfig(;kernel = :gaussian),
	RbfConfig(;kernel = :multiquadric),
	RbfConfig(;kernel = :inv_multiquadric),
	RbfConfig(;kernel = :cubic),
	RbfConfig(;kernel = :gaussian, polynomial_degree = -1),
	RbfConfig(;kernel = :inv_multiquadric, polynomial_degree = -1),
	RbfConfig(;kernel = :cubic, polynomial_degree = -1),
	RbfConfig(;kernel = :gaussian, polynomial_degree = 0),
	RbfConfig(;kernel = :inv_multiquadric, polynomial_degree = 0),
	RbfConfig(;kernel = :cubic, polynomial_degree = 0),
	TaylorConfig(;degree=1),
	TaylorConfig(;degree=2),
	TaylorApproximateConfig(mode=:autodiff),
	TaylorApproximateConfig(mode=:fdm)
]

# TODO `TaylorCallbackConfig` cannot be tested this way

NVAR_RANGE = 1:3
NOUT_RANGE = 1:2
#%%

function perform_tests( cfg, n_vars, n_out, k, unconstrained, other_func_indices, x0 = Float64[] )
	f1 = scalar_test_funcs[k]

	lb = fill(-5.0, n_vars )
	ub = fill(6.0, n_vars )

	if unconstrained
		mop = Morbit.MixedMOP(n_vars)
	else
		mop = Morbit.MixedMOP(lb,ub)
	end

	test_func = isempty(other_func_indices) ? f1 : x -> [ f1(x); [ φ(x) for φ in scalar_test_funcs[other_func_indices] ] ]
	Morbit.add_vector_objective!( mop, test_func, cfg; n_out )

	if isempty(x0) x0 = Morbit._rand_box_point(lb, ub) end
	@show x0

	smop, iter_data, data_base, sc, ac = Morbit.initialize_data(mop, x0, Float32[] )

	_mod = sc.surrogates[1].model
	objf = Morbit.list_of_objectives(smop)[1]
	meta = sc.surrogates[1].meta

	# `Morbit.eval_objf` evaluates at a site from the original domain
	_eval_objf_original = x -> Morbit.eval_objf(objf, x)
	_eval_objf_scaled = x -> Morbit.eval_objf( objf, Morbit.unscale(x, smop) )

	# `eval_models`, `get_gradient` etc., on the other hand,
	# take input from [0,1]^n:
	_eval_model = x -> Morbit.eval_models( _mod, x )

	x0_s = Morbit.scale( x0, smop )

	@test x0_s ≈ Morbit.get_x( iter_data )
	@test _eval_objf_original( x0 ) ≈ _eval_objf_scaled( x0_s )

	@test isapprox(Morbit.get_gradient( _mod, x0_s, 1), ForwardDiff.gradient( x -> _eval_model(x)[1], x0_s ); rtol = 0.001 )
	@test isapprox(Morbit.get_jacobian( _mod, x0_s ), ForwardDiff.jacobian( _eval_model, x0_s ); rtol = 0.001 )

	for i = 1 : 10
		ξ0= Morbit._rand_box_point( lb, ub )
		ξ0_s = Morbit.scale( ξ0, smop )
		@test isapprox(Morbit.get_gradient( _mod, ξ0_s, 1), ForwardDiff.gradient( x -> _eval_model(x)[1], ξ0_s ); rtol = 0.001 )
		@test isapprox(Morbit.get_jacobian( _mod, ξ0_s ), ForwardDiff.jacobian( _eval_model, ξ0_s ); rtol = 0.001 )
	end
end#function

for cfg in configs
	for n_vars = NVAR_RANGE
		for n_out = NOUT_RANGE
			for k in eachindex(scalar_test_funcs)
				for unconstrained in [false,true]
					@info n_vars, n_out, k, unconstrained
					@info cfg

					other_func_indices = rand( collect(eachindex(scalar_test_funcs)), n_out - 1 )

#					for other_func_indices in combinations( eachindex(scalar_test_funcs), n_out - 1 )
						perform_tests(cfg, n_vars, n_out, k, unconstrained, other_func_indices )
#					end
				end#for
			end#for
		end#for
	end#for
end#for

end#module

using .ModelDerivatives


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

function int64_to_bytes(i64::Int64)::Vector{UInt8}
    vcat(
        UInt8(i64 & 0xFF),
        UInt8(i64 >> 8 & 0xFF),
        UInt8(i64 >> 16 & 0xFF),
        UInt8(i64 >> 24 & 0xFF),
        UInt8(i64 >> 32 & 0xFF),
        UInt8(i64 >> 40 & 0xFF),
        UInt8(i64 >> 48 & 0xFF),
        UInt8(i64 >> 56 & 0xFF),
    )
end

function uint32_to_bytes(ui32::UInt32)::Vector{UInt8}
    vcat(
        UInt8(ui32 & 0xFF),
        UInt8(ui32 >> 8 & 0xFF),
        UInt8(ui32 >> 16 & 0xFF),
        UInt8(ui32 >> 24 & 0xFF),
    )
end

function int32_to_bytes(i32::Int32)::Vector{UInt8}
    vcat(
        UInt8(i32 & 0xFF),
        UInt8(i32 >> 8 & 0xFF),
        UInt8(i32 >> 16 & 0xFF),
        UInt8(i32 >> 24 & 0xFF),
    )
end

function bint64_to_bytes(i64::Int64)::Vector{UInt8}
    vcat(
        UInt8(i64 >> 56 & 0xFF),
        UInt8(i64 >> 48 & 0xFF),
        UInt8(i64 >> 40 & 0xFF),
        UInt8(i64 >> 32 & 0xFF),
        UInt8(i64 >> 24 & 0xFF),
        UInt8(i64 >> 16 & 0xFF),
        UInt8(i64 >> 8 & 0xFF),
        UInt8(i64 & 0xFF),
    )
end

function bint32_to_bytes(i32::Int32)::Vector{UInt8}
    vcat(
        UInt8(i32 >> 24 & 0xFF),
        UInt8(i32 >> 16 & 0xFF),
        UInt8(i32 >> 8 & 0xFF),
        UInt8(i32 & 0xFF),
    )
end

function bint16_to_bytes(i16::Int16)::Vector{UInt8}
    vcat(
        UInt8(i16 >> 24 & 0xFF),
        UInt8(i16 >> 16 & 0xFF),
        UInt8(i16 >> 8 & 0xFF),
        UInt8(i16 & 0xFF),
    )
end

function bytes_to_uint128(b::Vector{UInt8})::UInt128
    reinterpret(UInt128, b)[1]
end

function bytes_to_int64(b::Vector{UInt8})::Int64
    reinterpret(Int64, b)[1]
end

function bytes_to_uint32(b::Vector{UInt8})::UInt32
    reinterpret(UInt32, b)[1]
end

function bytes_to_int32(b::Vector{UInt8})::Int32
    reinterpret(Int32, b)[1]
end

function bytes_to_uint16(b::Vector{UInt8})::UInt16
    reinterpret(UInt16, b)[1]
end

function bytes_to_int16(b::Vector{UInt8})::Int16
    reinterpret(Int16, b)[1]
end

function bytes_to_buint128(b::Vector{UInt8})::UInt128
    reinterpret(UInt128, reverse(b))[1]
end

function bytes_to_buint64(b::Vector{UInt8})::UInt64
    reinterpret(UInt64, reverse(b))[1]
end

function bytes_to_buint32(b::Vector{UInt8})::UInt32
    reinterpret(UInt32, reverse(b))[1]
end

function bytes_to_buint16(b::Vector{UInt8})::UInt16
    reinterpret(UInt16, reverse(b))[1]
end

function bytes_to_bint128(b::Vector{UInt8})::Int128
    reinterpret(Int128, reverse(b))[1]
end

function bytes_to_bint64(b::Vector{UInt8})::Int64
    reinterpret(Int64, reverse(b))[1]
end

function bytes_to_bint32(b::Vector{UInt8})::Int32
    reinterpret(Int32, reverse(b))[1]
end

function bytes_to_bint16(b::Vector{UInt8})::Int16
    reinterpret(Int16, reverse(b))[1]
end

function bytes_to_str(b::Vector{UInt8})::String
    String(b)
end

function xdr_bytes(bs::Vector{UInt8})::Vector{UInt8}
    # XDR encoding bytes
    n = length(bs)
    padding = 0
    if n % 4 != 0
        padding = 4 - n % 4
    end
    buf = zeros(UInt8, 4 + n + padding)
    buf[1] = UInt8(n >> 24 & 0xFF)
    buf[2] = UInt8(n >> 16 & 0xFF)
    buf[3] = UInt8(n >> 8 & 0xFF)
    buf[4] = UInt8(n & 0xFF)
    for i in 1:length(bs)
        buf[4 + i] = bs[i]
    end
    buf
end

function to_blr(i64::Int64)::Tuple{Vector{UInt8}, Vector{UInt8}}
    bint64_to_bytes(i64), UInt8[16, 0]
end

function to_blr(i32::Int32)::Tuple{Vector{UInt8}, Vector{UInt8}}
    bint32_to_bytes(i32), UInt8[8, 0]
end

function to_blr(f64::Float64)::Tuple{Vector{UInt8}, Vector{UInt8}}
    buf = IOBuffer()
    write(buf, f64)
    seek(buf, 0)
    v = read(buf)
    v, UInt8[27]
end

function to_blr(bytes::Vector{UInt8})::Tuple{Vector{UInt8}, Vector{UInt8}}
    nbytes = length(bytes)
    pad_length = ((4 - nbytes) & 3)
    padding = Vector{UInt8}([0, 0, 0])
    v = vcat(bytes, padding[1:pad_length])

    v, UInt8[14, UInt8(nbytes & 255), UInt8(nbytes >> 8)]
end

function _convert_date(year::Int, month::Int, day::Int)::Vector{UInt8}
    i = month + 9
    jy = year + div(i, 12) - 1
    jm = i % 12
    c = div(jy, 100)
    jy -= 100 * c
    j = (146097*c)/4 + (1461*jy)/4 + (153*jm+2)/5 + day - 678882

    bint32_to_bytes(Int32(j))
end

function _convert_time(hour::Int, minute::Int, second::Int, microsecond::Int, nanosecond::Int)::Vector{UInt8}
    v = (hour*3600+mintes*60+second)*10000 + div(microsecond, 100) + div(nanosecond, 100000)

    bint32_to_bytes(Int32(v))
end

function to_blr(d::Date)::Tuple{Vector{UInt8}, Vector{Uint8}}
    v = _convert_date(Dates.year(d), Dates.month(d), Dates.day(d))

    v, UInt8[12]
end

function to_blr(t::Time)::Tuple{Vector{UInt8}, Vector{UInt8}}
    v = _convert_time(Time.hour(t), Time.minute(t), Time.second(t), Time.microsecond(t), Time.nanosecond(t))

    v, UInt8[13]
end

function to_blr(dt::DateTime)::Tuple{Vector{UInt8}, Vector{UInt8}}
    v = vcat(
        _convet_date(DateTime.year(dt), DateTime.month(dt), DateTime.day(dt)),
        _convet_time(DateTime.hour(dt), DateTime.minute(dt), DateTime.second(dt), DateTime.microsecond(dt), DateTime.nanosecond(dt))
    )
    v, UInt8[35]
end

function to_blr(b::Bool)::Tuple{Vector{UInt8}, Vector{UInt8}}
    if b
        UInt8[1, 0, 0, 0], UInt8[23]
    else
        UInt8[0, 0, 0, 0], UInt8[23]
    end
end

function to_blr(Nothing)::Tuple{Vector{UInt8}, Vector{UInt8}}
    [], UInt8[14, 0, 0]
end

function calc_blr(xsqlda::Vector{XSQLVAR})::Vector{UInt8}
    # Calculate  BLR from XSQLVAR array.
    ln = length(xsqlda) * 2
    blr::Vector{UInt8} = UInt8[5, 2, 4, 0, UInt8(ln & 255), UInt8(ln >> 8)]

    for x in xsqlda
        sqlscale = x.sqlscale
        if sqlscale < 0
            sqlscale += 256
        end
        if x.sqltype == SQL_TYPE_VARYING
            append!(blr, UInt8(37))
            append!(blr, UInt8(x.sqllen & 255))
            append!(blr, UInt8(x.sqllen >> 8))
        elseif x.sqltype == SQL_TYPE_TEXT
            append!(blr, UInt8(14))
            append!(blr, UInt8(x.sqllen & 255))
            append!(blr, UInt8(x.sqllen >> 8))
        elseif x.sqltype == SQL_TYPE_LONG
            append!(blr, UInt8(8))
            append!(blr, UInt8(sqlscale))
        elseif x.sqltype == SQL_TYPE_SHORT
            append!(blr, UInt8(7))
            append!(blr, UInt8(sqlscale))
        elseif x.sqltype == SQL_TYPE_INT64
            append!(blr, UInt8(16))
            append!(blr, UInt8(sqlscale))
        elseif x.sqltype == SQL_TYPE_INT128
            append!(blr, UInt8(26))
            append!(blr, UInt8(sqlscale))
        elseif x.sqltype == SQL_TYPE_QUAD
            append!(blr, UInt8(9))
            append!(blr, UInt8(sqlscale))
        elseif x.sqltype == SQL_TYPE_DEC_FIXED  # OBSOLATED
            append!(blr, UInt8(26))
            append!(blr, UInt8(sqlscale))
        elseif x.sqltype == SQL_TYPE_DOUBLE
            append!(blr, UInt8(27))
        elseif x.sqltype == SQL_TYPE_FLOAT
            append!(blr, UInt8(10))
        elseif x.sqltype == SQL_TYPE_D_FLOAT
            append!(blr, UInt8(11))
        elseif x.sqltype == SQL_TYPE_DATE
            append!(blr, UInt8(12))
        elseif x.sqltype == SQL_TYPE_TIME
            append!(blr, UInt8(13))
        elseif x.sqltype == SQL_TYPE_TIMESTAMP
            append!(blr, UInt8(35))
        elseif x.sqltype == SQL_TYPE_BLOB
            append!(blr, UInt8(9))
            append!(blr, UInt8(0))
        elseif x.sqltype == SQL_TYPE_ARRAY
            append!(blr, UInt8(9))
            append!(blr, UInt8(0))
        elseif x.sqltype == SQL_TYPE_BOOLEAN
            append!(blr, UInt8(23))
        elseif x.sqltype == SQL_TYPE_DEC64
            append!(blr, UInt8(24))
        elseif x.sqltype == SQL_TYPE_DEC128
            append!(blr, UInt8(25))
        elseif x.sqltype == SQL_TYPE_TIME_TZ
            append!(blr, UInt8(28))
        elseif x.sqltype == SQL_TYPE_TIMESTAMP_TZ
            append!(blr, UInt8(29))
        end

        # [blr_short, 0]
        append!(blr, UInt8(7))
        append!(blr, UInt8(0))
    end

    # [blr_end, blr_eoc]
    append!(blr, UInt8(255))
    append!(blr, UInt8(76))

    blr
end


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

using dummyMindReader

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

#  parse shell arguments
shArgs = shArgParser(ARGS)

begin
  file = shArgs["file"]
  output = shArgs["output"]
end

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

mindReader(file, output)

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

max_stars_repo_name: DanielRivasMD/dummyMindReader

max_issues_repo_name: DanielRivasMD/dummyMindReader

max_forks_repo_name: DanielRivasMD/dummyMindReader


text:
using Altro
using TrajectoryOptimization
using Test
if !isdefined(Main,:TEST_TIME)
    TEST_TIME = true
end

# Double Integrator
solver = ALTROSolver(Problems.DoubleIntegrator()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 1
@test max_violation(solver) < 1e-6
@test iterations(solver) == 8 # 8
@test solver.stats.gradient[end] < 1e-9
@test status(solver) == Altro.SOLVE_SUCCEEDED

# Pendulum
solver = ALTROSolver(Problems.Pendulum()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 2
@test max_violation(solver) < 1e-6
@test iterations(solver) == 19 # 19
@test solver.stats.gradient[end] < 1e-3
@test status(solver) == Altro.SOLVE_SUCCEEDED

# Cartpole
solver = ALTROSolver(Problems.Cartpole()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 <  10
@test max_violation(solver) < 1e-6
@test iterations(solver) == 40 # 40
@test solver.stats.gradient[end] < 1e-2
@test status(solver) == Altro.SOLVE_SUCCEEDED

solver = ALTROSolver(Problems.Cartpole()..., projected_newton=false)
if !Sys.iswindows()
    @test b = benchmark_solve!(solver).allocs == 0
end

solver = ALTROSolver(Problems.Cartpole()..., projected_newton=false, static_bp=false)
if !Sys.iswindows()
    @test b = benchmark_solve!(solver).allocs == 0
end

# Acrobot
solver = ALTROSolver(Problems.Acrobot()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 10
@test max_violation(solver) < 1e-6
@test iterations(solver) == 50 # 50
@test solver.stats.gradient[end] < 1e-2
@test status(solver) == Altro.SOLVE_SUCCEEDED

# Parallel Park
solver = ALTROSolver(Problems.DubinsCar(:parallel_park)...)
b =  benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time /1e6 < 7
@test max_violation(solver) < 1e-6
@test iterations(solver) == 13 # 13
@test solver.stats.gradient[end] < 1e-3
@test status(solver) == Altro.SOLVE_SUCCEEDED

# Three Obstacles
solver = ALTROSolver(Problems.DubinsCar(:three_obstacles)...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time /1e6 < 6
@test max_violation(solver) < 1e-6
@test iterations(solver) == 20 # 20
@test solver.stats.gradient[end] < 1e-2
@test status(solver) == Altro.SOLVE_SUCCEEDED

solver = ALTROSolver(Problems.DubinsCar(:three_obstacles)..., projected_newton=false)
@test solver.opts.projected_newton == false
@test solver.stats.gradient[end] < 1e-1
if !Sys.iswindows()   # not sure why this fails on Windows?
    @test benchmark_solve!(solver).allocs == 0
    @test status(solver) == Altro.SOLVE_SUCCEEDED
end

# Escape
solver = ALTROSolver(Problems.DubinsCar(:escape)..., infeasible=true, R_inf=0.1)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 20
@test max_violation(solver) < 1e-6
@test iterations(solver) == 13 # 13
@test solver.stats.gradient[end] < 1e-3
@test status(solver) == Altro.SOLVE_SUCCEEDED

# Zig-zag
solver = ALTROSolver(Problems.Quadrotor(:zigzag)...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 60
@test max_violation(solver) < 1e-6
@test iterations(solver) == 15 # 16
@test solver.stats.gradient[end] < 2e-2
@test status(solver) == Altro.SOLVE_SUCCEEDED

solver = ALTROSolver(Problems.Quadrotor(:zigzag)..., projected_newton=false)
@test solver.stats.gradient[end] < 0.3
if !Sys.iswindows()   # not sure why this fails on Windows?
    @test benchmark_solve!(solver).allocs == 0
    @test status(solver) == Altro.SOLVE_SUCCEEDED
end

solver = ALTROSolver(Problems.Quadrotor(:zigzag)..., projected_newton=false,
    infeasible=true, static_bp=false, constraint_tolerance=1e-4)
b = benchmark_solve!(solver, samples=2, evals=2)
if !Sys.iswindows()
    @test b.allocs == 0
end
@test iterations(solver) == 20
@test status(solver) == Altro.SOLVE_SUCCEEDED

# Barrell Roll
solver = ALTROSolver(Problems.YakProblems()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 100
@test max_violation(solver) < 1e-6
@test iterations(solver) == 18 # 18
@test solver.stats.gradient[end] < 1e-3
@test status(solver) == Altro.SOLVE_SUCCEEDED


text:
module MaxHelpingHandBadelineSprite

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/BadelineSprite" BadelineSprite(x::Integer, y::Integer, left::Bool=false, floating::Bool=false)

const placements = Ahorn.PlacementDict(
    "Badeline Sprite (max480's Helping Hand)" => Ahorn.EntityPlacement(
        BadelineSprite
    )
)

function Ahorn.selection(entity::BadelineSprite)
    x, y = Ahorn.position(entity)

    scaleX = get(entity.data, "left", false) ? -1 : 1

    spriteName = "characters/badeline/idle00"
    sprite = Ahorn.getSprite(spriteName, "Gameplay")

    return Ahorn.getSpriteRectangle(spriteName, x, y, jx=0.5, jy=1.0, sx=scaleX)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BadelineSprite, room::Maple.Room)
    scaleX = get(entity.data, "left", false) ? -1 : 1
    spriteName = "characters/badeline/idle00"

    Ahorn.drawSprite(ctx, spriteName, 0, 0, jx=0.5, jy=1.0, sx=scaleX)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using JuliaConSamplePackage
using Test

@testset "JuliaConSamplePackage.jl" begin
    @test JuliaConSamplePackage.newton(sin, cos, 1.0) ≈ 0.0
    @test JuliaConSamplePackage.newton(cos, x -> -sin(x), 1.0) ≈ pi/2
    @test JuliaConSamplePackage.newton(x -> x^2 - 2x - 1, x -> 2x - 2, 3.0)  ≈ 2.414213562373095
    @test abs(JuliaConSamplePackage.newton(x -> exp(x) - cos(x), x -> exp(x) + sin(x), 3.0) - 0.0) <= 1e-14
    @test abs(JuliaConSamplePackage.newton(x -> x^5 - x - 1, x -> 5x^4 - 1, 1.0) - 1.1673039782614187) <= 1e-14

    @testset "autodiff" begin
        @test abs(JuliaConSamplePackage.newton(sin, 0.5) - 0.0) <= 100*eps(1.0)
        @test JuliaConSamplePackage.newton(cos, 1.0) ≈ pi/2
        @test JuliaConSamplePackage.newton(x -> x^2 - 2x - 1, 3.0)  ≈ 2.414213562373095
        @test abs(JuliaConSamplePackage.newton(x -> exp(x) - cos(x), 3.0) - 0.0) <= 1e-14
        @test abs(JuliaConSamplePackage.newton(x -> cos(x) - x^3, 3.0) - 0.8654740331016144) <= 1e-14
    end
end


text:
function system_model!(out, dx, x, (controls, sys), t)

    #Index Setup
    bus_size = get_bus_size(sys)
    bus_range = get_bus_range(sys)
    bus_vars_count = length(bus_range)
    injection_start = get_injection_pointer(sys)
    injection_count = 1
    branches_start = get_branches_pointer(sys)
    branches_count = 1

    #Network quantities
    V_r = @view x[1:bus_size]
    V_i = @view x[bus_size+1:bus_vars_count]
    Sbase = get_sys_base(sys)
    I_injections_r = zeros(bus_size)
    I_injections_i = zeros(bus_size)
    injection_ode = zeros(get_n_injection_states(sys))
    branches_ode = zeros(get_n_branches_states(sys))

    for d in sys.dyn_injections
        bus_n = get_bus_number(d)
        ix_range = range(injection_start, length=d.n_states)
        ode_range = range(injection_count, length=d.n_states)
        injection_count = injection_count + d.n_states
        injection_start = injection_start + d.n_states
        device_model!(x,
                      injection_ode,
                      view(V_r, bus_n),
                      view(V_i, bus_n),
                      view(I_injections_r, bus_n),
                      view(I_injections_i, bus_n),
                      ix_range,
                      ode_range,
                      controls,
                      d,
                      sys)
        out[ix_range] = injection_ode[ode_range] - dx[ix_range]
    end
    if !isnothing(sys.dyn_branch)
        for br in sys.dyn_branch
            arc = br.arc
            n_buses = get_bus_size(sys)
            from_bus_number = PSY.get_number(arc.from)
            to_bus_number = PSY.get_number(arc.to)
            ix_dx = [from_bus_number,
                    from_bus_number+n_buses,
                    to_bus_number,
                    to_bus_number+n_buses]
            ix_range = range(branches_start, length=br.n_states)
            ode_range = range(branches_count, length=br.n_states)
            branches_count = branches_count + br.n_states
            branch_model!(x,
                        dx,
                        branches_ode,
                        #Get Voltage data
                        view(V_r, from_bus_number),
                        view(V_i, from_bus_number),
                        view(V_r, to_bus_number),
                        view(V_i, to_bus_number),
                        #Get Current data
                        view(I_injections_r, from_bus_number),
                        view(I_injections_i, from_bus_number),
                        view(I_injections_r, to_bus_number),
                        view(I_injections_i, to_bus_number),
                        ix_range,
                        ix_dx,
                        ode_range,
                        controls,
                        br,
                        sys)

            out[ix_range] = branches_ode[ode_range] - dx[ix_range]
        end
    end

    for d in sys.injections
        bus_n = get_bus_number(d)

        device_model!( view(V_r, bus_n),
                       view(V_i, bus_n),
                       view(I_injections_r, bus_n),
                       view(I_injections_i, bus_n),
                       d,
                       sys)
    end

    out[bus_range] = kcl(sys.Ybus, V_r, V_i,
                               I_injections_r, I_injections_i)

end


text:
######################################################################
# sample_initializer.jl
######################################################################

"""
`si = SampleInitializer()`

Samples points from a distribution.
"""
mutable struct SampleInitializer <: Initializer end

function initialize(si::SampleInitializer, em::ErgodicManager, tm::TrajectoryManager)
	xd = Array(Vector{Float64}, tm.N+1)
	points = Array(Vector{Float64}, tm.N)
	xd[1] = deepcopy(tm.x0)
	ud = Array(Vector{Float64}, tm.N)
	bin_size = (em.bins, em.bins)

	# first sample N points from e.phi
	weights = Weights(vec(em.phi))
	for n = 1:tm.N
		xi, yi = ind2sub(bin_size, sample(weights))
		points[n] = [em.cell_size*(xi-.5), em.cell_size*(yi-.5)]
	end

	# find a short path heuristically
	#tsp_rand!(xd, points)
	tsp_nn!(xd, points)

	# compute controls
	ud = compute_controls(xd, tm.h)

	return xd, ud
end


text:
function logmvbeta(vals)
    r = 0.0
    for v in vals
        r += lgamma(v)
    end
    r -= lgamma(sum(vals))

    return r
end

function λ_objective(λ::Vector{Float64}, ∇λ::Vector{Float64},
        ν::Vector{Float64}, Ndivζ::Vector{Float64}, sumθ::Vector{Float64},
        μ::Vector{Float64}, invΣ::Matrix{Float64})

    diff = λ .- μ
    Eeη = exp.(λ .+ 0.5ν)

    if length(∇λ) > 0
        ∇λ .= -invΣ * diff .+ sumθ .- Ndivζ .* Eeη
    end

    return -0.5 * (diff' * invΣ * diff)[1] + sum(λ .* sumθ) - sum(Ndivζ .* Eeη)
end

function ν_objective(ν::Vector{Float64}, ∇ν::Vector{Float64},
        λ::Vector{Float64}, Ndivζ::Vector{Float64}, μ::Vector{Float64},
        invΣ::Matrix{Float64})

    Eeη = exp.(λ .+ 0.5ν)

    if length(∇ν) > 0
        ∇ν .= -0.5diag(invΣ) .- (Ndivζ / 2) .* Eeη .+ (1 ./ (2ν))
    end

    return -0.5 * tr(diagm(0 => ν) * invΣ) - sum(Ndivζ .* Eeη) + sum(log.(ν)) / 2
end

function α_objective(α::Vector{Float64}, ∇α::Vector{Float64},
        sum_Elnϕ::Float64, K::Int, V::Int)

    if length(∇α) > 0
        ∇α[1] = K * V * (digamma(V * α[1]) - digamma(α[1])) + sum_Elnϕ
    end

    return K * (lgamma(V * α[1]) - V * lgamma(α[1])) + α[1] * sum_Elnϕ
end

function check_convergence(metric::Vector{Vector{Float64}}; tol=1e-4)
    reldiff = maximum(abs.(metric[end - 1] .- metric[end]) ./ abs.(metric[end]))
    return reldiff < tol
end

function check_convergence(metric::Vector{Float64}; tol=1e-4)
    reldiff = maximum(abs(metric[end - 1] - metric[end]) / abs(metric[end]))
    return reldiff < tol
end


text:
using PyPlot

function visPreprocessing(KS::KineticSolver,file::String)

	rho = readdlm(file,',')
	#println("Final mass = $(sum(rho[:]))." )
	#nicewrite(string(KS.outputfolderprefix,"data/rhofinal.txt"),rho)

	if KS.ictype == "linesource"
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = 0.0; maxVal = 0.5;
	#	title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude), rot=$(KS.rotationmagnitude)", fontsize=20)
		pcolormesh(rho, cmap="RdBu_r",vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho.png"))

		rhoana = readdlm("../vis/exactLineSource.txt",',')
		#fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		plot(rhoana[:,1],rhoana[:,2])
		rhong  = rho[3:end-2,3:end-2] # remove ghosts
		cuthori = rhong[Int64(KS.ny/2),:]
		x = range(-1.5,stop=1.5,length=KS.nx-4)
		plot(x,cuthori)
		cutverti = rhong[:,Int64(KS.nx/2)]
		plot(x,cutverti)
		cutdiag = [rhong[i,i] for i=1:KS.nx-4]
		plot(sqrt(2)*x,cutdiag)
		cutdiag = [rhong[KS.ny-4+1-i,i] for i=1:KS.nx-4]
		plot(sqrt(2)*x,cutdiag)
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), convolution=$(KS.convolutionmagnitude) ", fontsize=30)
		PyPlot.savefig(string(KS.outputfolderprefix,"rhocut.png"))

	elseif KS.ictype == "checkerboard"
		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(log10.(abs.(rho)), cmap="RdBu_r",vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho_log10.png"))

	elseif KS.ictype == "beam"

		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(log10.(abs.(rho[3:end-2,3:end-2])), cmap="RdBu_r")#,vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho_log10.png"))

		fig, ax = subplots(figsize=(10.5, 8), dpi=100)
		minVal = -7.0; maxVal = 0.0;
		title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)", fontsize=20)
		pcolormesh(((rho[3:end-2,3:end-2])), cmap="RdBu_r")#,vmin=minVal, vmax=maxVal)
		colorbar()
		PyPlot.savefig(string(KS.outputfolderprefix,"rho.png"))

	end
end

max_stars_repo_name: camminady/SN

max_issues_repo_name: camminady/SN

max_forks_repo_name: camminady/SN


text:
@inline function two_sum(a::T, b::T) where T<:AbstractFloat
    ab =  a + b
    B  = ab - a
    ε  = (a - (ab - B)) + (b - B)
    return ab, ε
end

@inline function two_diff(a::T, b::T) where T<:AbstractFloat
    ab =  a - b
    B  = ab - a
    ε  = (a - (ab - B)) + (b - B)
    return ab, ε
end

"""
    fast_two_sum(x, y)::(a, b)

{ ab  =  a + b,  ε  = (a - ab) + b }

__unchecked precondition__: abs(a) >= abs(b).
"""
@inline function fast_two_sum(a::T, b::T) where T<:AbstractFloat
     ab =  a + b
     B  = ab - a
     ε  =  b - B
     return ab, ε
end

"""
    fast_two_diff(x, y)::(a, b)

{ ab  =  a - b,  ε  = (a - ab) + b }

__unchecked precondition__: abs(a) >= abs(b).
"""
@inline function fast_two_diff(a::T, b::T) where T<:AbstractFloat
     ab =  a - b
     B  = ab - a
     ε  =  B - b
     return ab, ε
end

@inline function two_prod(a::T, b::T) where T<:AbstractFloat
    ab =  a * b
    ε  = fma(a, b, -ab)
    return ab, ε
end

#=
   Separate an IEEEFloat (fp) into two parts, (hi, lo), assuring
      `fp === hi + lo`, `hi === fp - lo`, `lo === fp - hi`, and
      `ulp(abs(hi)) > ufp(abs(lo))`.

      ulp, ufp are acronyms for u[nit in the] {l[ast], f[irst]} p[lace]

      ulp(fp) = eps(fp)/2
      ufp(fp) = ldexp(0.5, exponent(fp)

   ref
       A floating-point technique for extending the available precision.
       T. J. Dekker.
       Numer. Math., 18:224–242, 1971
=#

splitter(::Type{Float64}) =
    Float64(one(Int64) + one(Int64) << (cld(precision(Float64),2)-1))
splitter(::Type{Float32}) =
    Float32(one(Int32) + one(Int32) << (cld(precision(Float32),2)-1))
splitter(::Type{Float16}) =
    Float16(one(Int16) + one(Int16) << (cld(precision(Float16),2)-1))

splitmax(::Type{Float64}) = realmax(Float64) / splitter(Float64)
splitmax(::Type{Float32}) = realmax(Float32) / splitter(Float32)
splitmax(::Type{Float16}) = realmax(Float16) / splitter(Float16)

# Veldkamp splitting of a floating point value
@inline function splitting(x::T) where T<:AbstractFloat
    (!isfinite(x) || abs(x) > splitmax(T)) && throw(OverflowError("$x overflows"))
    z   = x * splitter(T)
    zₕᵢ = z - (z - x)
    zₗₒ = x - zₕᵢ
    return zₕᵢ, zₗₒ
end

# Rump splitting
#    In extractscalar, a floating-point number fp is split relative to p2,
#    a fixed power of 2.
function extractscalar(fp, p2=2^27)
    hi = (p2 + fp) - p2
    lo = fp - hi
    return hi, lo
end

#=
    splitfp(fp) =?= extractscalar(2^(?+exponent(fp)), fp)

=#

for (U,F) in ((:UInt64, :Float64), (:UInt32, :Float32), (:UInt16, :Float16))
  @eval begin
    @inline function ufp(x::$F)
        u = reinterpret($U, x)
        u = (u >> (precision($F)-1)) << (precision($F)-1)
        return reinterpret($F, u)
    end
  end
end

const Float64ulp = inv(ldexp(1.0, precision(Float64)))
const Float32ulp = inv(ldexp(1.0, precision(Float32)))
const Float16ulp = inv(ldexp(1.0, precision(Float16)))

@inline ulp(x::Float64) = ufp(x) * Float64ulp
@inline ulp(x::Float32) = ufp(x) * Float32ulp
@inline ulp(x::Float16) = ufp(x) * Float16ulp

const Float64eps = inv(ldexp(1.0, precision(Float64)-1))
const Float32eps = inv(ldexp(1.0, precision(Float32)-1))
const Float16eps = inv(ldexp(1.0, precision(Float16)-1))

@inline epsi(x::Float64) = ufp(x) * Float64eps
@inline epsi(x::Float32) = ufp(x) * Float32eps
@inline epsi(x::Float16) = ufp(x) * Float16eps


text:
using Random
using LinearAlgebra
using Plots
using Statistics

#Step 1: Initialization
n=1000 #Impact limit
N=50 #Number of nodes
t=100000
P=60
I = [rand(1:n) for i in 1:N] #Node impact vector
Imean=[sum(I)/N]
display(I)
popularity=[]
TopConnections=0

HighestImpact=[(findmax(I))[1]]

C = bitrand((N,N)) #Network matrix
display(C)
Connections=[sum(C)]

#plotlyjs()
function StateReport(stage)

    plot0=scatter(I,[sum(C[i,:]) for i in 1:N],legend=false)
    xlims!(0,1000)
    xlabel!("Impact")
    ylims!(0,N)
    ylabel!("Connectivity")
    title!(string("State",string(stage)))
    display(plot0)
    png(string(stage))
end

#Energy
E=[]
Enode = zeros(N,t)
function computeE()
    #println("Computing E")
    for i in 1:N
        B=0
        for k in 1:N
            if C[i,k] == 1
                B += I[k]
            end
        end

        Enode[i] = I[i]*B
    end
    push!(E,sum(Enode))
end

#Step 2: Network dynamics
p=0.2 #Event probability
#Remove lowest impact connection
#gaseste rand((1,N))
#loop and find lowest
function Step(i)
    minVal=n
    min=0
    candidates=[]
    #println("Node ",i)
    Execute1=rand(1:100)
    Execute2=rand(1:100)

    for j in 1:N
        if (C[i,j]==1) && (minVal>I[j])
            minVal=I[j]
            min=j
        end
    end

    if min!=0 && Execute1>P
        C[i,min]=0
        #C.data[min,i]=0
        #println("   Removed ",min," of impact ",minVal)
    end

    for j in 1:N
        if I[j]>minVal && C[i,j]==0 #! Se pot deconecta de ei insisi?
            push!(candidates,j)
        end
    end

    if length(candidates)>0 && Execute2>P
        j=rand(1:length(candidates))
        C[i,j]=1
        #C.data[j,i]=1
        #println("   Added ",j," of impact ",I[j])
    end

    #Hai sa nu se deconexcteze de ei insisi
    ##O fi bine aici??
    for i in 1:N
        C[i,i]=1
    end

    I[(findmin(I))[2]]=rand(1:n)
    #println("I[",absMin[2],"] cu val ",absMin[1]," inlocuit de ",I[absMin[2]])

    #Ar putea fi și în afara funcției
    #absMax=[(findmax(I))[2]]
    global TopConnections=0
    global Connections
    global HighestImpact
    for i in 1:N
        if I[i]==(findmax(I))[1]
            #push!(absMax,(findmax(I))[2])
            global TopConnections+=sum(C[:,i])
        end
    end
    #global Connections=sum(C[(findmax(I))[2],:])
    global popularity
    push!(popularity,TopConnections#=/length(absMax)=#)

    computeE()
    push!(Imean,sum(I)/N)
    push!(Connections,sum(C))
    push!(HighestImpact,(findmax(I))[1])

end

computeE()
E0=pop!(E)
StateReport(0)

for i in 1:t
    Step(rand((1:N)))
    if mod(i,1000)==0
        println("t=",i)
    end
end

plot1=plot(1:t,E/E0,legend=false)
xlabel!("Time")
ylabel!("E/E0")
xlims!(51000,52000)
#ylims!(2,2.5)
title!("Energy")
display(plot1)
png("plot1")

plot2=plot(1:(t+1),Imean,legend=false)
xlabel!("Time")
ylabel!("Mean impact")
xlims!(51000,51500)
ylims!(950,1000)
title!("Mean Impact over Time")
display(plot2)
png("plot2b")

plot3=plot(1:t,popularity,legend=false)
xlabel!("Time")
ylabel!("# connections for max. impact agents")
title!("Top impact popularity")
display(plot3)
png("plot3")

plot4=plot(1:(t+1),Connections,legend=false)
xlabel!("Time")
ylabel!("Total connections")
title!("Total connections")
display(plot4)
png("plot4")

plot5=plot(1:(t+1),HighestImpact,legend=false)
xlabel!("Time")
ylabel!("Highest impact")
title!("Highest impact over time")
display(plot5)
png("plot5")

#########################
###Degree distribution###
#########################
DegreeArray=[sum(C[i,:]) for i in 1:N]
plot6=histogram(DegreeArray, nbins = 30)
title!("Degree distribution")
display(plot6)
png("plot6")

####################
###Hurst Exponent###
####################
#=

#1. Calculate the mean
m=mean(popularity)
#2. Create a mean adjusted series
Y=[popularity[i]-m for i in 1:t]
#display(Y)
#3. Calculate the cumulative deviate series Z
Z=[sum(Y[1:i]) for i in 1:t]
#display(Z)
#4. Compute the range R
R=[maximum(Z[1:i])-minimum(Z[1:i]) for i in 1:t]
#display(R)
#5. Compute the standard deviation S
S=[std(popularity[1:i]) for i in 1:t]
#display(S)
#6. plotting log[R(n)/S(n)] as a function of log n
RescaledRange=[]
for i in 1:t
    global RescaledRange
    push!(RescaledRange, [R[i]/S[i]])
end
#display(RescaledRange[1000:t])
plot7=scatter(1000:t,RescaledRange[1000:t],legend=false)
xaxis!(log)
yaxis!(log)
xlabel!("log n")
ylabel!("log[R(n)/S(n)]")
title!("Hurst exponent")
display(plot7)
png("plot7")=#

max_stars_repo_path: Scripts/Prototype 1/prototype1.jl

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_path: Scripts/Prototype 1/prototype1.jl

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_path: Scripts/Prototype 1/prototype1.jl

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
using Images, FileIO
using ImageView

include("settings.jl")
include("dna.jl")
include("renderer.jl")
include("fitness.jl")

orig_img_path = "./monalisa.jpg" # ARGS[1]
orig_img = load(orig_img_path)

(MAX_WIDTH, MAX_HEIGHT) = size(orig_img)

current_pic = Pic()

generation = 0
selected = 0
error_level = Inf

while true
    global orig_img
    global MAX_WIDTH
    global MAX_HEIGHT

    global current_pic

    global generation
    global selected
    global error_level

    new_pic = deepcopy(current_pic)
    mutate!(new_pic)

    if new_pic.dirty
        generation += 1

        img = render(new_pic, 1, MAX_WIDTH, MAX_HEIGHT, generation)
        new_error_level = fitness(img, orig_img)

        if new_error_level <= error_level
            selected += 1

            current_pic = new_pic
            error_level = new_error_level

            if selected % 100 == 0
                @info "generation = $generation"
                @info "selected = $selected"
                @info "error_level = $error_level"
            end
        end
    end
end

max_stars_repo_name: llrt/evolisa-jl

max_issues_repo_name: llrt/evolisa-jl

max_forks_repo_name: llrt/evolisa-jl


text:
abstract type AbstractIPSolver end
abstract type AbstractIPOptions end

function interior_point_options(ip_type::Symbol)
    if ip_type == :interior_point
        return Symbol("InteriorPointOptions")
    elseif ip_type == :mehrotra
        return Symbol("InteriorPointOptions")
    else
        error("Unknown ip_type.")
    end
end

mutable struct ResidualMethods
    r!
    rz!
    rθ!
end

# residual
function r!(r, z, θ, κ)
    @warn "residual not defined"
    error()
    nothing
end

# residual Jacobian wrt z
function rz!(rz, z, θ)
    @warn "residual Jacobian wrt z not defined"
    nothing
end

# residual Jacobian wrt θ
function rθ!(rθ, z, θ)
    @warn "residual Jacobian wrt θ not defined"
    nothing
end

# optimization spaces
abstract type Space end

# Euclidean
struct Euclidean <: Space
    n::Int
end

function candidate_point!(z̄::Vector{T}, ::Euclidean, z::Vector{T}, Δ::Vector{T}, α::T) where T
    z̄ .= z - α .* Δ
end

function update_point!(z::Vector{T}, ::Space, z̄::Vector{T}) where T
    z .= z̄
end

function mapping!(δz, s::Euclidean, δzs, z) # TODO: make allocation free
    δz .= δzs
end

# interior-point solver options
@with_kw mutable struct InteriorPointOptions{T} <: AbstractIPOptions
    r_tol::T = 1.0e-5
    κ_tol::T = 1.0e-5
    ls_scale::T = 0.5
    max_iter::Int = 100
    max_ls::Int = 3
    max_time::T = 1e5
    diff_sol::Bool = false
    reg::Bool = false
    ϵ_min = 0.05 # ∈ [0.005, 0.25]
        # smaller -> faster
        # larger  -> slower, more robust
    κ_reg = 1e-3 # bilinear constraint violation level at which regularization is triggered [1e-3, 1e-4]
    γ_reg = 1e-1 # regularization scaling parameters ∈ [0, 0.1]:
        # 0   -> faster & ill-conditioned
        # 0.1 -> slower & better-conditioned
        # simulation choose γ_reg = 0.1
        # MPC choose γ_reg = 0.0
    solver::Symbol = :lu_solver
    undercut::T = 5.0 # the solver will aim at reaching κ_vio = κ_tol / undercut
        # simulation choose undercut = Inf
        # MPC choose undercut = 5.0
    verbose::Bool = false
    warn::Bool = false
end

mutable struct InteriorPoint{T,R,RZ,Rθ} <: AbstractIPSolver
    s::Space
    oss::OptimizationSpace
    methods::ResidualMethods
    z::Vector{T}               # current point
    r::R                       # residual
    rz::RZ                     # residual Jacobian wrt z
    rθ::Rθ                     # residual Jacobian wrt θ
    Δ::Vector{T}               # search direction
    δz::Matrix{T}              # solution gradients (this is always dense)
    δzs::Matrix{T}             # solution gradients (in optimization space; δz = δzs for Euclidean)
    θ::Vector{T}               # problem data
    solver::LinearSolver
    reg_val::T
    iterations::Int
    opts::InteriorPointOptions{T}
    κ
end

function interior_point(z, θ;
        s = Euclidean(length(z)),
        oss = OptimizationSpace12(),
        r! = r!, rz! = rz!, rθ! = rθ!,
        r  = zeros(oss.nΔ),
        rz = spzeros(oss.nΔ, oss.nΔ),
        rθ = spzeros(oss.nΔ, length(θ)),
        opts::InteriorPointOptions = InteriorPointOptions()) where T

    rz!(rz, z, θ) # compute Jacobian for pre-factorization
    num_data = length(θ)

    TYPE = typeof.([z[1], r, rz, rθ])
    InteriorPoint{TYPE...}(
        s,
        oss,
        ResidualMethods(r!, rz!, rθ!),
        z,
        r,
        rz,
        rθ,
        zeros(oss.nΔ),
        zeros(oss.nz, num_data),
        zeros(oss.nΔ, num_data),
        θ,
        eval(opts.solver)(rz),
        0.0,
        0,
        opts,
        zeros(1),
        )
end

# interior point solver
function interior_point_solve!(ip::InteriorPoint{T,R,RZ,Rθ}) where {T,R,RZ,Rθ}

    # space
    s = ip.s
    oss = ip.oss

    # options
    opts = ip.opts
    r_tol = opts.r_tol
    κ_tol = opts.κ_tol
    ls_scale = opts.ls_scale
    max_iter = opts.max_iter
    max_time = opts.max_time
    max_ls = opts.max_ls
    diff_sol = opts.diff_sol
    ϵ_min = opts.ϵ_min
    κ_reg = opts.κ_reg
    γ_reg = opts.γ_reg
    reg = opts.reg
    verbose = opts.verbose
    warn = opts.warn

    # unpack pre-allocated data
    z = ip.z
    r = ip.r
    rz = ip.rz
    Δ = ip.Δ
    θ = ip.θ

    # indices
    ortz = oss.ortz
    ortΔ = oss.ortΔ
    socz = oss.socz
    socΔ = oss.socΔ
    bil = oss.bil
    ortr = oss.ortr
    socr = oss.socr

    # Initialization
    solver = ip.solver
    ip.iterations = 0
    no_progress = 0

    # compute residual, residual Jacobian
    ip.methods.r!(r, z, θ, 0.0)
    # TODO need least squares init /!|
    # least_squares!(ip, z, θ, r, rz) # seems to be harmful for performance (failure and iteration count)
    z .= initial_state!(z, ortz, socz) # decrease failure rate for linearized case

    ip.methods.r!(r, z, θ, 0.0)

    κ_vio = bilinear_violation(ip, r)
    r_vio = residual_violation(ip, r)
    elapsed_time = 0.0

    for j = 1:max_iter
        elapsed_time >= max_time && break
        elapsed_time += @elapsed begin
            # check for converged residual
            if (r_vio < r_tol) && (κ_vio < κ_tol)
                break
            end
            ip.iterations += 1

            # Compute regularization level
            κ_vio = bilinear_violation(ip, r)
            ip.reg_val = κ_vio < κ_reg ? κ_vio * γ_reg : 0.0

            # compute residual Jacobian
            rz!(ip, rz, z, θ, reg = ip.reg_val) # this is not adapted to the second order cone

            # compute step
            linear_solve!(solver, Δ, rz, r, reg = ip.reg_val)

            α_ort = ort_step_length(z, Δ, ortz, ortΔ; τ = 1.0)
            α_soc = soc_step_length(z, Δ, socz, socΔ; τ = 1.0, verbose = false)
            α = min(α_ort, α_soc)
            μ, σ = general_centering(z, Δ, ortz, ortΔ, socz, socΔ, α)
            αaff = α

            # Compute corrector residual
            ip.methods.r!(r, z, θ, max(σ * μ, κ_tol/opts.undercut)) # here we set κ = σ*μ, Δ = Δaff
            general_correction_term!(r, Δ, ortr, socr, ortΔ, socΔ)

            # Compute corrector search direction
            linear_solve!(solver, Δ, rz, r, reg = ip.reg_val, fact = false)
            τ = max(0.95, 1 - max(r_vio, κ_vio)^2)

            α_ort = ort_step_length(z, Δ, ortz, ortΔ; τ = τ)
            α_soc = soc_step_length(z, Δ, socz, socΔ; τ = min(τ, 0.99), verbose = false)
            α = min(α_ort, α_soc)

            # reduce norm of residual
            candidate_point!(z, s, z, Δ, α)
            κ_vio_cand = 0.0
            r_vio_cand = 0.0
            for i = 1:max_ls
                ip.methods.r!(r, z, θ, 0.0)
                κ_vio_cand = bilinear_violation(ip, r)
                r_vio_cand = residual_violation(ip, r)
                if (r_vio_cand <= r_vio) || (κ_vio_cand <= κ_vio)
                    break
                end
                verbose && println("linesearch $i")
                # backtracking
                candidate_point!(z, s, z, Δ, -α * ls_scale^i)
            end

            κ_vio = κ_vio_cand
            r_vio = r_vio_cand
            # nc = Int(oss.ny/4)
            # nb = Int(oss.ny/2)
            verbose && println("iter:", j,
                "  r: ", scn(norm(r, Inf)),
                "  r_vio: ", scn(r_vio),
                "  κ_vio: ", scn(κ_vio),
                # "  κ: ", scn(norm(r[oss.bil[1:nc]], Inf)),
                # "  κ: ", scn(norm(r[oss.bil[nc .+ (1:nb)]], Inf)),
                # "  κ: ", scn(norm(r[oss.bil[nc + nb .+ (1:nc)]], Inf)),
                "  Δ: ", scn(norm(Δ)),
                "  α: ", scn(norm(α)))

            # verbose && println(
            #     "in:", j,
            #     "   αaff:", scn(αaff, digits = 0),
            #     "   α:", scn(α, digits = 0),
            #     "   μσ:", scn(μ*σ, digits = 0),
            #     "   κ_vio:", scn(κ_vio, digits = 0),
            #     "   r_vio:", scn(r_vio, digits = 0),
            #     )
        end
    end
    if (r_vio < r_tol) && (κ_vio < κ_tol)
        # differentiate solution
        ########################################################################
        # @warn "changed"
        if false && R <: RLin
        # if R <: RLin
            # regularize solution so that k_vio == k_tol for all bilinear constraints
            ip.methods.r!(r, z, θ, 0.0)
            for i in eachindex(oss.ortz[1])
                if z[oss.ortz[1][i]] <= z[oss.ortz[2][i]]
                    z[oss.ortz[1][i]] *= κ_tol / r.rbil[i]
                else
                    z[oss.ortz[2][i]] *= κ_tol / r.rbil[i]
                end
            end
        else
            # # regularize solution so that k_vio == k_tol for all bilinear constraints
            # ip.methods.r!(r, z, θ, 0.0)
            # for i in eachindex(oss.ortz[1])
            #     if z[oss.ortz[1][i]] <= z[oss.ortz[2][i]]
            #         z[oss.ortz[1][i]] *= κ_tol / r[oss.ortr][i]
            #     else
            #         z[oss.ortz[2][i]] *= κ_tol / r[oss.ortr][i]
            #     end
            # end
            # ip.methods.r!(r, z, θ, 0.0)
            # @show scn.(r[oss.ortr])
        end
        ########################################################################
        diff_sol && differentiate_solution!(ip, reg = max(ip.reg_val, κ_tol * γ_reg))
        return true
    else
        return false
    end
end

function rz!(ip::AbstractIPSolver, rz::AbstractMatrix{T}, z::AbstractVector{T},
        θ::AbstractVector{T}; reg = 0.0) where {T}
    z_reg = deepcopy(z)
    ortz = ip.oss.ortz
    socz = ip.oss.socz
    for i in eachindex(ortz) # primal-dual
        z_reg[ortz[i]] = max.(z[ortz[i]], reg)
    end
    ip.methods.rz!(rz, z_reg, θ)
    return nothing
end

function rθ!(ip::AbstractIPSolver, rθ::AbstractMatrix{T}, z::AbstractVector{T},
        θ::AbstractVector{T}) where {T}
    ip.methods.rθ!(rθ, z, θ)
    return nothing
end

function general_correction_term!(r::AbstractVector{T}, Δ::AbstractVector{T},
    ortr::Vector{Int}, socr::Vector{Int},
    ortΔ::Vector{Vector{Int}}, socΔ::Vector{Vector{Vector{Int}}}) where {T}
    # @warn "define residual order"
    r[ortr] .+= Δ[ortΔ[1]] .* Δ[ortΔ[2]]
    r[socr] .+= vcat(
        [second_order_cone_product(
            Δ[socΔ[2][i]],
            Δ[socΔ[1][i]],
        ) for i in eachindex(socΔ[1])]...)
    return nothing
end

function least_squares!(ip::AbstractIPSolver, z::AbstractVector{T}, θ::AbstractVector{T},
        r::AbstractVector{T}, rz::AbstractMatrix{T}) where {T}
    # doing nothing gives the best result if z_t is correctly initialized with z_t-1 in th simulator
        dyn = ip.oss.dyn
        rst = ip.oss.rst
        A = rz[[dyn; rst], :]
        z .+= A' * ((A * A') \ r[[dyn; rst]])
    return nothing
end

function initial_state!(z::AbstractVector{T}, ortz::Vector{Vector{Int}},
        socz::Vector{Vector{Vector{Int}}}; ϵ::T=1e-20) where {T}

    # Split between primals and duals
    socz_p = socz[1]
    socz_d = socz[2]
    ortz_p = ortz[1]
    ortz_d = ortz[2]

    # ineq
    y1 = z[ortz_p]
    y2 = z[ortz_d]
    δy1 = max(-1.5 * minimum(y1), 0)
    δy2 = max(-1.5 * minimum(y2), 0)

    y1h = y1 .+ δy1
    y2h = y2 .+ δy2
    δhy1 = 0.5 * y1h'*y2h / (sum(y2h) + ϵ)
    δhy2 = 0.5 * y1h'*y2h / (sum(y1h) + ϵ)

    y10 = y1h .+ δhy1
    y20 = y2h .+ δhy2
    z[ortz_p] .= y10
    z[ortz_d] .= y20

    # soc
    for i in eachindex(socz_p)
        e = [1; zeros(length(socz_p[i]) - 1)] # identity element
        y1 = z[socz_p[i]]
        y2 = z[socz_d[i]]
        δy1 = max(-1.5 * (y1[1] - norm(y1[2:end])), 0)
        δy2 = max(-1.5 * (y2[1] - norm(y2[2:end])), 0)

        y1h = y1 + δy1 * e
        y2h = y2 + δy2 * e
        δhy1 = 0.5 * y1h'*y2h / ((y2h[1] + norm(y2h[2,end])) + ϵ)
        δhy2 = 0.5 * y1h'*y2h / ((y1h[1] + norm(y1h[2,end])) + ϵ)

        y10 = y1h + δhy1 * e
        y20 = y2h + δhy2 * e
        z[socz_p[i]] .= y10
        z[socz_d[i]] .= y20
    end
    return z
end

function interior_point_solve!(ip::AbstractIPSolver, z::AbstractVector{T}, θ::AbstractVector{T}) where T
    ip.z .= z
    ip.θ .= θ
    interior_point_solve!(ip)
end

function differentiate_solution!(ip::AbstractIPSolver; reg = 0.0)
    s = ip.s
    z = ip.z
    θ = ip.θ
    rz = ip.rz
    rθ = ip.rθ
    δz = ip.δz
    δzs = ip.δzs

    rz!(ip, rz, z, θ, reg = reg)
    rθ!(ip, rθ, z, θ)

    linear_solve!(ip.solver, δzs, rz, rθ, reg = reg)
    @inbounds @views @. δzs .*= -1.0
    mapping!(δz, s, δzs, z)
    return nothing
end

################################################################################
# New methods
################################################################################

function residual_violation(ip::AbstractIPSolver, r::AbstractVector{T}) where {T}
    dyn = ip.oss.dyn
    rst = ip.oss.rst
    max(norm(r[dyn], Inf), norm(r[rst], Inf))
end

function general_centering(z::AbstractVector{T}, Δaff::AbstractVector{T},
        ortz::Vector{Vector{Int}}, ortΔ::Vector{Vector{Int}},
        socz::Vector{Vector{Vector{Int}}}, socΔ::Vector{Vector{Vector{Int}}}, αaff::T) where {T}
        # See Section 5.1.3 in CVXOPT
        # μ only depends on the dot products (no cone product)
        # The CVXOPT linear and quadratic cone program solvers

    n = length(ortz[1]) + sum(length.(socΔ[1]))
    # ineq
    μ = z[ortz[1]]' * z[ortz[2]]
    μaff = (z[ortz[1]] - αaff * Δaff[ortΔ[1]])' * (z[ortz[2]] - αaff * Δaff[ortΔ[2]])
    # soc
    for i in eachindex(socz[1])
        μ += z[socz[1][i]]' * z[socz[2][i]]
        μaff += (z[socz[1][i]] - αaff * Δaff[socΔ[1][i]])' * (z[socz[2][i]] - αaff * Δaff[socΔ[2][i]])
    end
    μ /= n
    μaff /= n
	σ = clamp(μaff / μ, 0.0, 1.0)^3
	return μ, σ
end

function bilinear_violation(ip::AbstractIPSolver, r::AbstractVector{T}) where {T}
    bil = ip.oss.bil
    return norm(r[bil], Inf)
end

function soc_value(u::AbstractVector{T}) where {T}
    u0 = u[1]
    u1 = u[2:end]
    return (u0^2 - u1' * u1)
end

function soc_step_length(λ::AbstractVector{T}, Δ::AbstractVector{T};
        τ::T = 0.99, ϵ::T = 1e-14, verbose::Bool = false) where {T}
    # check Section 8.2 CVXOPT
    # The CVXOPT linear and quadratic cone program solvers

    # Adding to slack ϵ to make sure that we never get out of the cone
    λ0 = λ[1] #- ϵ
    λ_λ = max(λ0^2 - λ[2:end]' * λ[2:end], 1e-25)
    verbose && println(
        "    vλ:", scn(soc_value(λ), digits = 0, exp_digits = 2),
        "    vλ+Δ:", scn(soc_value(λ+Δ), digits = 0, exp_digits = 2),
        "    λ_λ: ", scn(λ_λ, digits = 0, exp_digits = 2),
        "    λ:", scn.(λ, digits = 0, exp_digits = 2),
        "    Δ:", scn.(Δ, digits = 0, exp_digits = 2),
        )
    if λ_λ < 0.0
        @show λ_λ
        @warn "should always be positive"
        # error("should always be positive")
    end
    λ_λ += ϵ
    λ_Δ = λ0 * Δ[1] - λ[2:end]' * Δ[2:end] + ϵ

    ρs = λ_Δ / λ_λ
    ρv = Δ[2:end] / sqrt(λ_λ)
    ρv -= (λ_Δ / sqrt(λ_λ) + Δ[1]) / (λ0 / sqrt(λ_λ) + 1) * λ[2:end] / λ_λ
    # we make sre that the inverse always exists with ϵ,
    # if norm(ρv) - ρs) is negative (Δ is pushing towards a more positive cone)
        # the computation is ignored and we get the maximum value for α = 1.0
    # else we have α = τ / norm(ρv) - ρs)
    # we add ϵ to the denumerator to ensure strict positivity and avoid 1e-16 errors.
    α = 1.0
    if norm(ρv) - ρs > 0.0
        α = min(α, τ / (norm(ρv) - ρs))
    end
    verbose && println(
        "     α:", scn(α, digits = 0, exp_digits = 2))
    verbose && cone_plot(λ, Δ, show_α = true)
    return α
end

function soc_step_length(z::AbstractVector{T}, Δ::AbstractVector{T},
		socz::Vector{Vector{Vector{Int}}}, socΔ::Vector{Vector{Vector{Int}}};
        τ::T=0.99, verbose::Bool = false) where {T}
        # We need to make this much more efficient (allocation free)
    α = 1.0
    for i in eachindex(socz) # primal-dual
        for j in eachindex(socz[i]) # number of cones
            # we need -Δ here because we will taking the step x - α Δ
            α = min(α, soc_step_length(z[socz[i][j]], -Δ[socΔ[i][j]], τ = τ, verbose = verbose))
        end
    end
    return α
end

function ort_step_length(z::AbstractVector{T}, Δ::AbstractVector{T},
		ortz::Vector{Vector{Int}}, ortΔ::Vector{Vector{Int}};
        τ::T=0.9995) where {T}
        # We need to make this much more efficient (allocation free)
    α = 1.0
    for i in eachindex(ortz) # primal-dual
        for j in eachindex(ortz[i])
            k = ortz[i][j] # z
            ks = ortΔ[i][j] # Δz
            if Δ[ks] > 0.0
                α = min(α, τ * z[k] / Δ[ks])
            end
        end
    end
    return α
end

#650 lines before cleanup


text:
using DiffEqSensitivity
using OrdinaryDiffEq, Calculus, Test

function f(du,u,p,t)
  du[1] = -p[1]
  du[2] = p[2]
end

cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(println("Stopped.");integrator.p[2]=0.0))
p = [2.0, 1.0]

function test_f(p)
  prob = ODEProblem(f,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
  solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,callback=cb,save_everystep=false)[end]
end
p = [2.0, 1.0]
findiff = Calculus.finite_difference_jacobian(test_f,p)
findiff

using ForwardDiff
ad = ForwardDiff.jacobian(test_f,p)
ad

@test ad ≈ findiff


text:
module CliqueInstanceReader

function readInstance(filename, cut_of_no_nodes)
    isfile(filename) || error("cannot find data-file: $(filename)")
    ff = open(filename)

    no_nodes::Int32=0
    no_edges::Int32=0
    while( !eof(ff) )
        line = readline(ff)
        if line[1]=='c'
            #println(li, " c: $(line)")
        elseif line[1]=='p'
            #println("p $(line)")
            data = split(line)
            no_nodes= parse(Int32,data[3])
            no_edges= parse(Int32,data[4])
            break
        elseif line[1]=='e'
            li+=1
            exit(0)
        else
            println("wrong format: >$(line)<")
            exit(0)'
        end
    end

    if cut_of_no_nodes<0
       cut_of_no_nodes=no_nodes
    end
    println("datafile: ",filename," nodes: ", no_nodes, " edges: ", no_edges)
    println("cut_of_no_nodes: ",cut_of_no_nodes)
    if no_nodes<cut_of_no_nodes
        cut_of_no_nodes=no_nodes
    end
    edges=zeros(Int8,cut_of_no_nodes,cut_of_no_nodes)
    while( !eof(ff) )
        line = readline(ff)
        if line[1]=='c'
            #println(li, " c: $(line)")
        elseif line[1]=='p'
            #println(li, "p $(line)")
        elseif line[1]=='e'
            data = split(line)
            la_node= parse(Int32,data[2])
            sm_node= parse(Int32,data[3])
            if la_node<=cut_of_no_nodes  && sm_node<=cut_of_no_nodes
                edges[la_node,sm_node]=1
                edges[sm_node,la_node]=1
            end
        else
            println("wrong format: >$(line)<")
            exit(0)'
        end
    end

    return edges
end

export readInstance

end

max_stars_repo_name: henrikfoss/CliqueInstanceReaderPkg

max_issues_repo_name: henrikfoss/CliqueInstanceReaderPkg

max_forks_repo_name: henrikfoss/CliqueInstanceReaderPkg


text:
abstract type BuildPhase end

abstract type Existing <: BuildPhase end    #Existing projects.
abstract type Planned <: BuildPhase end     #Projects under construction.
abstract type Queue <: BuildPhase end       #Projects in queue.
abstract type Option <: BuildPhase end      #Projects which can be built.
abstract type Retired <: BuildPhase end     #Retired projects.


text:
# the better parts of this derive from ChangePrecision.jl

module SaferIntTypes

export @saferintegers

import Random
using Random: AbstractRNG

using MacroTools

using ..SaferIntegers

## Note: code in this module must be very careful with math functions,
#        because we've defined module-specific versions of very basic
#        functions like + and *.   Call Base.:+ etcetera if needed.

const HWInt = Union{Bool,Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128}

const SFInt = Union{SafeInt, SafeInt8, SafeInt16, SafeInt32, SafeInt64, SafeInt128, SafeUInt, SafeUInt8, SafeUInt16, SafeUInt32, SafeUInt64, SafeUInt128}

const randfuncs = (:rand,) # random-number generators
const matfuncs = (:ones, :zeros) # functions to construct arrays
const binaryfuncs = (:*, :+, :-, :^) # binary functions on irrationals that make Float64

const changefuncs = Set([randfuncs..., matfuncs...,
                         binaryfuncs..., :include])

changetype(x) = x

changetype(x::SafeInt128) = x
changetype(x::SafeInt64) = x
changetype(x::SafeInt32) = x
changetype(x::SafeInt16) = x
changetype(x::SafeInt8) = x

changetype(x::SafeUInt128) = x
changetype(x::SafeUInt64) = x
changetype(x::SafeUInt32) = x
changetype(x::SafeUInt16) = x
changetype(x::SafeUInt8) = x

changetype(x::Int128) = SafeInt128(x)
changetype(x::Int64) = SafeInt64(x)
changetype(x::Int32) = SafeInt32(x)
changetype(x::Int16) = SafeInt16(x)
changetype(x::Int8) = SafeInt8(x)

changetype(x::UInt128) = SafeUInt128(x)
changetype(x::UInt64) = SafeUInt64(x)
changetype(x::UInt32) = SafeUInt32(x)
changetype(x::UInt16) = SafeUInt16(x)
changetype(x::UInt8) = SafeUInt8(x)

changetype(x::String) = x

function changetype(ex::Expr)
    if isa(ex, String)
        return ex
    elseif isa(ex, Integer)
        return Int === Int64 ? SafeInt64(ex) : SafeInt32(ex)
    elseif (Meta.isexpr(ex, :call, 2) && isa(ex.args[2], Integer))
        return changetype_of_cast_int(ex)
    elseif (Meta.isexpr(ex, :(=), 2) && isa(ex.args[2], Integer) && isa(ex.args[1], Symbol))
        ex.args[2] = Int === Int64 ? SafeInt64(ex.args[2]) : SafeInt32(ex.args[2])
        return ex
    elseif (Meta.isexpr(ex, :(=), 2) && isa(ex.args[2], Integer))
        return changetype_of_cast_int(ex)
   # elseif (Meta.isexpr(ex, :call, 2) && isa(ex.args[2].args[2], Integer))
    #    return changetype_of_assigned_cast_int(ex)
    #elseif (Meta.isexpr(ex, :(=), 2) && isa(ex.args[2].args[2], Integer))
     #   return changetype_of_assigned_cast_int(ex)
    else
        return changetypes(ex)
    end
end

function changetypes(ex::Expr)
    if Meta.isexpr(ex, :call, 3) && ex.args[1] == :^ && ex.args[3] isa Int
        # mimic Julia 0.6/0.7's lowering to literal_pow
        return Expr(:call, ChangeType.literal_pow, :^, changetype(ex.args[2]), Val{ex.args[3]}())
    elseif Meta.isexpr(ex, :call, 2) && ex.args[1] == :include
        return :($include(@__MODULE__, $(ex.args[2])))
    elseif Meta.isexpr(ex, :call) && ex.args[1] in changefuncs
        return Expr(:call, Core.eval(SaferIntTypes, ex.args[1]), changetype.(ex.args[2:end])...)
    elseif Meta.isexpr(ex, :., 2) && ex.args[1] in changefuncs && Meta.isexpr(ex.args[2], :tuple)
        return Expr(:., Core.eval(SaferIntTypes, ex.args[1]), Expr(:tuple, changetype.(ex.args[2].args)...))
    elseif Meta.isexpr(ex, :call, 3) && ex.args[1] == :^ && ex.args[3] isa Int
    else
        return Expr(ex.head, changetype.(ex.args)...)
    end
end

function changetype_of_cast_int(ex::Expr)
   if @capture(ex, Int64(_))
      ex.args[1] = :SafeInt64
   elseif @capture(ex, Int32(_))
      ex.args[1] = :SafeInt32
   elseif @capture(ex, Int16(_))
      ex.args[1] = :SafeInt16
   elseif @capture(ex, Int8(_))
      ex.args[1] = :SafeInt8
   elseif @capture(ex, Int128(_))
      ex.args[1] = :SafeInt128
   elseif @capture(ex, UInt64(_))
      ex.args[1] = :SafeUInt64
   elseif @capture(ex, UInt32(_))
      ex.args[1] = :SafeUInt32
   elseif @capture(ex, UInt16(_))
      ex.args[1] = :SafeUInt16
   elseif @capture(ex, UInt8(_))
      ex.args[1] = :SafeUInt8
   elseif @capture(ex, UInt128(_))
      ex.args[1] = :SafeUInt128
   end
   return ex
end

function changetype_of_assigned_cast_int(ex::Expr)
   if @capture(ex, _=Int64(_))
      ex.args[2].args[1] = :SafeInt64
   elseif @capture(ex, _=Int32(_))
      ex.args[2].args[1] = :SafeInt32
   elseif @capture(ex, _=Int16(_))
      ex.args[2].args[1] = :SafeInt16
   elseif @capture(ex, _=Int8(_))
      ex.args[2].args[1] = :SafeInt8
   elseif @capture(ex, _=Int128(_))
      ex.args[2].args[1] = :SafeInt128
   elseif @capture(ex, _=UInt64(_))
      ex.args[2].args[1] = :SafeUInt64
   elseif @capture(ex, _=UInt32(_))
      ex.args[2].args[1] = :SafeUInt32
   elseif @capture(ex, _=UInt16(_))
      ex.args[2].args[1] = :SafeUInt16
   elseif @capture(ex, _=UInt8(_))
      ex.args[2].args[1] = :SafeUInt8
   elseif @capture(ex, _=UInt128(_))
      ex.args[2].args[1] = :SafeUInt128
   end
   return ex
end

# calls to include(f) are changed to include(T, mod, f) so that
# @changetype can apply recursively to included files.
function include(mod, filename::AbstractString)
    # use the undocumented parse_input_line function so that we preserve
    # the filename and line-number information.
    s = string("begin; ", read(filename, String), "\nend\n")
    expr = Base.parse_input_line(s, filename=filename)
    Core.eval(mod, changetype(expr))
end

for F in binaryfuncs
#    @eval $F(T, args...) = Base.$F(changetype.(T, args...)...,)
    @eval Base.$F(::Type{S}, ::Type{I}, n) where {S<:SFInt, I<:HWInt} = $F(S, n)
end

for F in randfuncs
    @eval Base.$F(::Type{S}, ::Type{I}, n) where {S<:SFInt, I<:HWInt} = $F(S, n)
end

for F in matfuncs
    @eval Base.$F(::Type{S}, ::Type{I}, n) where {S<:SFInt, I<:HWInt} = $F(S, n)
end

macro saferintegers(x)
   esc(changetype(x))
end

end # SaferIntTypes


text:
# unpacking and packing for optimization toolboxes.
function packbcW(b, c::Vector{T}, W) where T <: Real
    D = length(b[1])
    M = length(W)
    @assert length(b) == length(c) == M

    v = Vector{T}(undef,M+M*D+M*div(D*(D+1),2))

    st = 1
    fin = M
    v[st:fin] = c

    for m = 1:M
        v[fin+(m-1)*D+1:fin+m*D] = b[m]
    end
    st = fin+1
    fin = fin+M*D

    for m = 1:M
        L = cholesky(W[m]).L
        chol_c, chol_a = packCholeskyType(L)
        st = fin+1
        fin = fin + D
        v[st:fin] = chol_c

        st = fin+1
        fin = fin + div(D*(D+1),2)-D
        v[st:fin] = chol_a
    end

    return v
end

# forces c to be positive. Forces W to be posdef.
function unpackbcW(v::Vector{T}, M::Int, D::Int) where T <: Real

    st = 1
    fin = M
    c = v[st:fin]
    c = abs.(c) # force positive.

    b = collect( v[fin+(m-1)*D+1:fin+m*D] for m = 1:M )
    st = fin+1
    fin = fin+M*D

    W = Vector{Matrix{T}}(undef,M)
    for m = 1:M
        st = fin+1
        fin = fin + D
        chol_c = v[st:fin]
        chol_c = abs.(chol_c) # force positive.

        st = fin+1
        fin = fin + div(D*(D+1),2)-D
        chol_a = v[st:fin]

        L = parseCholeskyType(chol_c, chol_a)
        W[m] = L*L'
    end

    return b, c, W
end

# modifies L such that the lower triangular off-diagonal entries are replaced by a.
# The diagonal entries of L are unmodified.
function filloffdiagonalentriessym!(L::Matrix{T}, a::Vector{T})::Nothing where T <: Real
    N = size(L,1)
    @assert size(L,2) == size(L,1)

    p = 1
    for j = 1:N
        for i = j+1:N
            L[i,j] = a[p]
            p += 1
        end
    end
    @assert p-1 == round(Int, N*(N+1)/2-N) # debug.

    return nothing
end

function parseCholeskyType(c::Vector{T}, a::Vector{T})::Matrix{T} where T <: Real
    L = convert(Matrix{T}, LinearAlgebra.Diagonal(c))

    filloffdiagonalentriessym!(L, a)

    return L
end

# c is the diagonal, should be all positives.
function packCholeskyType(L)
    c = collect( L[i,i] for i = 1:size(L,1) )

    n = length(c)
    a = Vector{Float64}(undef, round(Int, n*(n+1)/2-n))

    p = 1
    for j = 1:n
        for i = j+1:n
            a[p] = L[i,j]
            p += 1
        end
    end
    @assert p-1 == round(Int, n*(n+1)/2-n) # debug.

    return c, a
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
# %% imports

include("../../Modules/save_param.jl")
include("critic_sim.jl")
using YAML

# %% get the results

results = load_param("Julia/Cartpole/Criticv2/results002.jld")
function results_to_genes(results::Dict{String,Any}, indiv::Int64)
    genes_hyp = Float64[]
    hyp = String[]
    dic_keys = collect(keys(results))
    filter!(x-> x!="weights", dic_keys)
    filter!(x-> x!="fit", dic_keys)
    for key in dic_keys
        push!(hyp, key)
        push!(genes_hyp, results[key][indiv])
    end
    genes_alpha = map(x->borne(x, 0.0, 10.0),results["weights"][indiv])
    matalpha = reshape(genes_alpha, (4,4))
    return matalpha, genes_hyp, hyp
end

# %% create the networks

indiv = 20

#actor
filepath = "Julia/CartPole/LIF_base/results/results_EI888.jld"
n_actor, _ = init_teacher_student(filepath, 19)

#critic default
arch_critic = [8,8,2]
params = YAML.load(open("Julia/Cartpole/Criticv2/cfg_default_critic.yml"))
n_critic = Network(arch_critic, params)

#tune the critic with results
matalpha, genes_hyp,hyp = results_to_genes(results, indiv)
assert_genes!(n_critic,genes_hyp, hyp)

td_errors = Float64[]

for epoch in 1:20
    @show epoch
    td_error, td_errors = play_episode(n_actor, n_critic, matalpha)
    @show maximum(td_errors)
    display(plot(td_errors))
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP