text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGBalanceLawDiscretizations
using Printf
using LinearAlgebra
using Logging
using CLIMA.MPIStateArrays

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

const ArrayTypes = (CLIMA.array_type(),)

const _nauxstate = 5
const _a_r, _a_θ, _a_ϕ, _a_f, _a_fdown = 1:_nauxstate
@inline function auxiliary_state_initialization!(aux, x, y, z, Rinner, Router)
  @inbounds begin
    r = hypot(x, y, z)
    θ = atan(y , x)
    ϕ = asin(z / r)
    aux[_a_r] = r
    aux[_a_θ] = θ
    aux[_a_ϕ] = ϕ

    # Exact integral
    a = 1 + sin(θ)^2 + sin(ϕ)^2
    aux[_a_f] = exp(-a * r^2) - exp(-a * Rinner^2)
    aux[_a_fdown] = exp(-a * Router^2) - exp(-a * r^2)
  end
end

@inline function integral_knl(val, Q, aux)
  @inbounds begin
    r, θ, ϕ = aux[_a_r], aux[_a_θ], aux[_a_ϕ]
    a = 1 + sin(θ)^2 + sin(ϕ)^2
    val[1] = -2r * a * exp(-a * r^2)
  end
end

using Test
function run(mpicomm, topl, ArrayType, N, FT, Rinner, Router)
  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = Topologies.cubedshellwarp,
                                         )

  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = 0,
                           flux! = (x...) -> (),
                           numerical_flux! = (x...) -> (),
                           auxiliary_state_length = _nauxstate,
                           auxiliary_state_initialization! = (x...) ->
                           auxiliary_state_initialization!(x..., Rinner,
                                                           Router),
                           numerical_boundary_flux! = (x...) -> ())

  Q = MPIStateArray(spacedisc)
  exact_aux = copy(spacedisc.auxstate)

  DGBalanceLawDiscretizations.indefinite_stack_integral!(spacedisc,
                                                         integral_knl, Q,
                                                         _a_f)
  DGBalanceLawDiscretizations.reverse_indefinite_stack_integral!(spacedisc,
                                                                 _a_fdown, _a_f)

  euclidean_distance(exact_aux, spacedisc.auxstate)
end

let
  CLIMA.init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  polynomialorder = 4

  base_Nhorz = 4
  base_Nvert = 2
  Rinner = 1//2
  Router = 1

  polynomialorder = 4

  expected_result = [6.228615762850257e-7
                     9.671308320438864e-9
                     1.5102832678375277e-10
                     2.359860999112363e-12]

  lvls = integration_testing ? length(expected_result) : 1

  @testset "$(@__FILE__)" for ArrayType in ArrayTypes
    for FT in (Float64,) #Float32)
      err = zeros(FT, lvls)
      for l = 1:lvls
        @info (ArrayType, FT, "sphere", l)
        Nhorz = 2^(l-1) * base_Nhorz
        Nvert = 2^(l-1) * base_Nvert
        Rrange = range(FT(Rinner); length=Nvert+1, stop=Router)
        topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
        err[l] = run(mpicomm, topl, ArrayType, polynomialorder, FT,
                     FT(Rinner), FT(Router))
        @test expected_result[l] ≈ err[l]
      end
      if integration_testing
        @info begin
          msg = ""
          for l = 1:lvls-1
            rate = log2(err[l]) - log2(err[l+1])
            msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
          end
          msg
        end
      end
    end
  end
end

nothing

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
#! /usr/bin/env julia

# Copyright 2020 John Hanley.
#
# 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.

struct PPM

    function PPM(size_px)
        println("P3")
        println(size_px, " ", size_px)
        println("255")
    end
end
    function plot(grey_value)
        v = grey_value
        println(v, " ", v, " ", v)
    end

function _get_px_res()
    return parse(UInt32, ENV["MSET_PX_RESOLUTION"])
end

function mandelbrot_set(xc, yc, sz, px_resolution=_get_px_res())
    PPM(px_resolution)
    step = (2 * sz) / (px_resolution - 1)

    for j in 0 : px_resolution-1
        for i in 0 : px_resolution-1
            x0 = xc - sz + step * i
            y0 = yc - sz + step * j
            plot(_cycles_to_escape(x0, y0))
        end
    end
end

function _cycles_to_escape(x0, y0, max_iter=255)
    x, y, i = 0.0, 0.0, 0
    while x * x + y * y <= 4 && i < max_iter
        x, y = x * x - y * y + x0, 2 * x * y + y0
        i += 1
    end
    return i
end

xc = parse(Float64, ARGS[1])
yc = parse(Float64, ARGS[2])
sz = parse(Float64, ARGS[3])
t0 = time()

mandelbrot_set(xc, yc, sz)

elapsed = round(time() - t0, digits=3)
println(stderr, elapsed, " seconds")

ext: jl


text:
using ProgressiveHedging
const PH = ProgressiveHedging

using Test

using Distributed

using Ipopt
using JuMP
using MathOptInterface
const MOI = MathOptInterface

using Suppressor
using TimerOutputs

macro optional_using(pkg)
    quote
        try
            using $pkg
        catch e
            if typeof(e) != ArgumentError
                rethrow(e)
            end
        end
    end
end

@optional_using(Xpress)

include("common.jl")

@testset "ProgressiveHedging" begin
    @testset "Scenario Tree" begin
        include("test_tree.jl")
    end
    @testset "Sanity Checks" begin
        include("test_sanity.jl")
    end
    @testset "Utils" begin
        include("test_utils.jl")
    end
    @testset "Subproblem" begin
        include("test_subproblem.jl")
    end
    @testset "Workers" begin
        include("test_worker.jl")
    end
    @testset "Penalty Parameters" begin
        include("test_penalty.jl")
    end
    @testset "Callbacks" begin
        include("test_callback.jl")
    end
    @testset "Setup" begin
        include("test_setup.jl")
    end
    @testset "Algorithm" begin
        include("test_algorithm.jl")
    end
    @testset "Distributed" begin
        include("test_distributed.jl")
    end
    @testset "Post Processing" begin
        include("test_post.jl")
    end
end


text:
function test_crls()
  ⪅(x,y) = (x ≈ y) || (x < y)
  crls_tol = 1.0e-5;

  for npower = 1 : 4
    (b, A, D, HY, HZ, Acond, rnorm) = test(40, 40, 4, npower, 0);  # No regularization.

    (x, stats) = crls(A, b);
    resid = norm(A' * (A*x - b)) / norm(b)
    @printf("CRLS: Relative residual: %8.1e\n", resid);
    @test(resid <= crls_tol);
    @test(stats.solved);

    λ = 1.0e-3;
    (x, stats) = crls(A, b, λ=λ);
    resid = norm(A' * (A*x - b) + λ * x) / norm(b)
    @printf("CRLS: Relative residual: %8.1e\n", resid);
    @test(resid <= crls_tol);
    @test(stats.solved);
  end

  # Test with preconditioning.
  A = rand(10, 6); b = rand(10);
  M = InverseLBFGSOperator(10, 4);
  for _ = 1 : 6
    s = rand(10);
    y = rand(10);
    push!(M, s, y);
  end

  (x, stats) = crls(A, b, M=M);
  resid = norm(A' * M * (A * x - b)) / sqrt(dot(b, M * b));
  @printf("CRLS: Preconditioned residual: %8.1e\n", resid);
  @test resid <= crls_tol;

  # test trust-region constraint
  (x, stats) = crls(A, b)

  radius = 0.75 * norm(x)
  (x, stats) = crls(A, b, radius=radius)
  @test(stats.solved)
  @test(abs(radius - norm(x)) <= crls_tol * radius)

  opA = LinearOperator(A)
  (xop, statsop) = crls(opA, b, radius=radius)
  @test(abs(radius - norm(xop)) <= crls_tol * radius)

  # Code coverage.
  (b, A, D, HY, HZ, Acond, rnorm) = test(40, 40, 4, 3, 0);
  (x, stats) = crls(Matrix(A), b);
  (x, stats) = crls(sparse(Matrix(A)), b);
  show(stats);

  # Test b == 0
  (x, stats) = crls(A, zeros(size(A,1)))
  @test x == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = [I; rand(1:10, 2, 3)]
  b = A * ones(Int, 3)
  (x, stats) = crls(A, b)
  @test stats.solved

  # Test A positive semi-definite
  radius = 10.
  m,n = 10,7
  U = qr(rand(m,m)).Q
  V = qr(rand(n,n)).Q
  S = [diagm(0 => [0, 1.0e-6, 1, 4, 20, 15, 1.0e5]) ; zeros(3,7)]
  A = U * S * V
  p = V[:,1]; b = A'\p;
  (x, stats) = crls(A, b, radius=radius)
  @test stats.solved
  @test stats.status == "zero-curvature encountered"
  @test norm(x) ⪅ radius
end

test_crls()


text:
using Zygote
zygotefun2(a,b) = sigmoid_fast(a + b) + LoopVectorization.relu(a * b) * tanh(a - b)
zygotefun3(a,b,c) = log(tanh_fast(a*b) + 1 + exp(b - c))
@testset "Zygote" begin
    @show @__LINE__

    for T ∈ (Float32,Float64)
        x = randn(T, 217); y = randn(T, 217); z = randn(T, 217);
        w = randn(T, 217)
        # Test 1-arg and `tanh_fast` vs `tanh`
        gref1 = gradient(x -> dot(w,  map(tanh,      x)), x);
        glv_1 = gradient(x -> dot(w, vmap(tanh_fast, x)), x);

        @test  only(gref1) ≈ only(glv_1)

        # Test 2 arguments
        # Test 3 arguments
        gref2 = gradient(xyz -> dot(w,  map(zygotefun2, xyz...)), (x, y));
        g_lv2 = gradient(xyz -> dot(w, vmap(zygotefun2, xyz...)), (x, y));

        @test all(map(≈, only(gref2), only(g_lv2)))

        gref3 = gradient(xyz -> dot(w,  map(zygotefun3, xyz...)), (x, y, z));
        g_lv3 = gradient(xyz -> dot(w, vmap(zygotefun3, xyz...)), (x, y, z));

        @test all(map(≈, only(gref3), only(g_lv3)))
    end
end


text:
using Documenter, DiskJockey

makedocs(
    modules = [DiskJockey],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == true, analytics = "UA-5472810-7"),
    sitename = "DiskJockey.jl Documentation",
    authors = "Ian Czekala and contributors.",
    pages = ["Home" => "index.md",
    "Introduction" => "dynamical_mass_intro.md",
    "Installation" => "installation.md",
    "Conventions" => "conventions.md",
    "Models" => "models.md",
    "Formats" => "formats.md",
    "RADMC Setup" => "RADMC3D_setup.md",
    "Cookbook" => "cookbook.md",
    "Priors" => "priors.md",
    "API" => "api.md",
    "Changelog" => "changelog.md"],
)

deploydocs(
    repo = "github.com/iancze/DiskJockey.git",
)


text:
# SPDX-License-Identifier: X11
# 2020-11-15 (Live)
# Retype (4pt, 7pt)

function main()
  T = parse(Int, readline())
  for i ∈ 1:T
    print("Case #$i: ")
    solve()
  end
end

function solve()
  N,K,S = parse.(Int, (split ∘ readline)())
  println(K - 1 + min(N + 1, (K - S) + (N + 1) - S))
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
"""
    vsubkron!(G, L, R)

Calculates `G[i, j] = G[i, j] - (∑ₘ L[i, m] R[j, m])` where `m` can be omitted,
meaning L and R can be Vectors. Used for local updates.
"""
function vsubkron!(G::Matrix{Float64}, L::Vector{Float64}, R::Vector{Float64})
    @turbo for k in eachindex(L), l in eachindex(R)
        G[k, l] -= L[k] * R[l]
    end
end

function vsubkron!(G::Matrix{Float64}, L::Matrix{Float64}, R::Matrix{Float64})
    @turbo for k in axes(L, 1), l in axes(R, 1), m in axes(L, 2)
        G[k, l] -= L[k, m] * R[l, m]
    end
end

function vsubkron!(G::CMat64, L::CVec64, R::CVec64)
    @turbo for k in eachindex(L), l in eachindex(R)
        G.re[k, l] -= L.re[k] * R.re[l]
    end
    @turbo for k in eachindex(L), l in eachindex(R)
        G.re[k, l] += L.im[k] * R.im[l]
    end
    @turbo for k in eachindex(L), l in eachindex(R)
        G.im[k, l] -= L.im[k] * R.re[l]
    end
    @turbo for k in eachindex(L), l in eachindex(R)
        G.im[k, l] -= L.re[k] * R.im[l]
    end
end

function vsubkron!(G::BlockDiagonal, L::Tuple, R::Tuple)
    @inbounds for b in eachindex(G.blocks)
        vsubkron!(G.blocks[b], L[b], R[b])
    end
end

"""
    vsub!(trg, ::UniformScaling, src, i, N)

Calculates `trg[m, n] = I[m, (i:N:end)[n]] - src[m, (i:N:end)[n]]` where
`n = eachindex(i:N:end)`.
"""
function vsub!(trg::FVec64, ::UniformScaling, src::FMat64, i::Int, N::Int)
    @turbo for j in eachindex(trg)
        trg[j] = - src[j, i]
    end
    @inbounds trg[i] += 1.0
    nothing
end

function vsub!(trg::CVec64, ::UniformScaling, src::CMat64, i::Int, N::Int)
    @turbo for j in eachindex(trg.re)
        trg.re[j] = - src.re[j, i]
    end
    @inbounds trg.re[i] += 1.0
    @turbo for j in eachindex(trg.im)
        trg.im[j] = - src.im[j, i]
    end
    nothing
end

function vsub!(trg::FMat64, ::UniformScaling, src::FMat64, i::Int, N::Int)
    @inbounds for (k, l) in enumerate(i:N:size(src, 2))
        @turbo for j in axes(trg, 1)
            trg[j, k] = - src[j, l]
        end
        trg[l, k] += 1.0
    end
    nothing
end

function vsub!(trg::Tuple, ::UniformScaling, src::BlockDiagonal, i::Int, N::Int)
    @inbounds for b in eachindex(trg)
        vsub!(trg[b], I, src.blocks[b], i, N)
    end
    nothing
end

"""
    vmul!(trg, M, src, i, N)

Calculates `trg[m, n] = ∑ₖ M[n, k] * src[(i:N:end)[k], m]` where
`n, k = eachindex(i:N:end)`.
"""
function vmul!(trg::FVec64, M::Float64, src::FMat64, i::Int, N::Int)
    @turbo for j in eachindex(trg)
        trg[j] = M * src[i, j]
    end
    nothing
end

function vmuladd!(trg::FVec64, M::Float64, src::FMat64, i::Int, N::Int)
    @turbo for j in eachindex(trg)
        trg[j] += M * src[i, j]
    end
    nothing
end

function vmul!(trg::CVec64, M::ComplexF64, src::CMat64, i::Int, N::Int)
    vmul!(trg.re, real(M), src.re, i, N)
    vmuladd!(trg.re, -imag(M), src.im, i, N)
    vmul!(trg.im, real(M), src.im, i, N)
    vmuladd!(trg.im, imag(M), src.re, i, N)
    nothing
end

function vmul!(trg::FMat64, M::FMat64, src::FMat64, i::Int, N::Int)
    r = i:N:size(src, 1)
    @inbounds for i in axes(trg, 1), j in eachindex(r)
        tij = 0.0
        for (k, l) in enumerate(r)
            tij += M[j, k] * src[l, i]
        end
        trg[i, j] = tij
    end
    nothing
end

function vmul!(trg::NTuple, M::Union{FVec64, CVec64}, src::BlockDiagonal, i::Int, N::Int)
    @inbounds for b in eachindex(trg)
        vmul!(trg[b], M[b], src.blocks[b], i, N)
    end
    nothing
end


text:
################################################################################
#  Copyright 2021, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

@testset "Steady State Process" begin

    @testset "2003 System" begin
        # Example from: Reliability of Safety-Critical System - Theory and
        # Applications (ch. 5)
        global λ = 5e-4u"1/hr"
        global μ = 1e-1u"1/hr"

        std²³ = include(joinpath(_MSS.BASE_DIR,"test/elements/2003_system.jl"))

        solve!(std²³, SteadyStateProcess())
        prb = [np[end] for np in _MSS.get_sprop(std²³, :prob)]
        sol = [ μ*(2*λ*λ + 2*λ*μ + μ*μ + λ*μ),
                3*λ*μ*(λ + μ),
                6*μ*λ*λ,
                6*λ*λ*λ] ./ (6*λ*λ*λ + 11*λ*λ*μ + 5*λ*μ*μ + μ*μ*μ + λ*μ*μ)

        @test isapprox(sum(prb), 1.0, rtol=1e-6)
        @test all(isapprox.(prb, sol, rtol=1e-6))
    end

    @testset "2003 System with Measurements Uncertainty" begin
        # Example from: Reliability of Safety-Critical System - Theory and
        # Applications (ch. 5)
        global λ = (5e-4 ± 1e-5)u"1/hr"
        global μ = (1e-1 ± 1e-2)u"1/hr"

        std²³ = include(joinpath(_MSS.BASE_DIR,"test/elements/2003_system.jl"))

        solve!(std²³, SteadyStateProcess())
        prb = [np[end] for np in _MSS.get_sprop(std²³, :prob)]
        sol = [ μ*(2*λ*λ + 2*λ*μ + μ*μ + λ*μ),
                3*λ*μ*(λ + μ),
                6*μ*λ*λ,
                6*λ*λ*λ] ./ (6*λ*λ*λ + 11*λ*λ*μ + 5*λ*μ*μ + μ*μ*μ + λ*μ*μ)

        @test isapprox(sum(prb), 1.0, rtol=1e-6)
        @test all(isapprox.(prb, sol, rtol=1e-6))

        @test isapprox(_MSM.uncertainty(sum(prb)), 0.0, atol=1e-6)
    end
end


text:
using Constants

time1 = true
time2 = true
tm = 0.0

if ion_viscosity
    dtm_ = 0.5 * dtm
    eps_visc_max = 1./3 * eps_visc_max #artificial viscosity
else
    dtm_ = dtm
end

function init_spacegrid(dr)
    r = zeros(Float64, nz)

    if(geom=="spherical")
        r[1] = L + rmin
        r0[:,nregions+1] = r0[:,nregions+1] + rmin #correct also r0 array
    else
        r[1] = 0.
    end

    for k = 2:nz
        r[k] = r[k-1] + dr #note that dr is negative
    end

    if(geom=="spherical")
        r[nz] = rmin
    end

    return r
end

function init_art_viscosity(geom)
    eps_visc = zeros(Float64, nz, neqi+1)

    if geom=="slab"
		for k = 1:nz
		   eps_visc[k,1:neqi+1] = eps_visc_max
		end
	else #spherical
		y1 = log10(eps_visc_max)
		y2 = log10(eps_visc_max / eps_compress)
		x1 = log10(r[nz])
		x2 = log10(r[1])
		b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
		a  = ( y1 - b ) / x1
		for k = 1:nz
			eps_visc[k,1:neqi+1] = r[k]^a * 10^b
		end
	end

    return eps_visc
end

function init_variables(hydro)

    T0 = zeros(Float64, nz, nspec)
    N0 = zeros(Float64, nz, nspec)
    V0 = zeros(Float64, nz, nspec)

	for j = 1:nspec
		for m = 1:nregions
			for k = 1:nz
				if (r[k] >= r0[j,m]) & (r[k] <= r0[j,m+1])
					T0[k,j] = temp0[j,m]
					N0[k,j] = den0[j,m]
					V0[k,j] = vel0[j,m]
				end
			end
		end
	end

    for k = 1:nz
    	for j = 1:nspec
			hydro.rho[k,j] = N0[k,j] * mi[j]
			hydro.u[k,j] = V0[k,j]
			hydro.p[k,j] = hydro.rho[k,j] / mi[j] * T0[k,j]
			hydro.T[k,j] = T0[k,j]
		end
    end

    return hydro
end

function init_predictor_corrector(hydro)
	hydro.U1D_p = copy(hydro.U1D)
	hydro.U1D_c = copy(hydro.U1D)
    return hydro
end

#----------------------------------------------------
function do_smoothing(hydro)

	kgrad = Array{Int32}(2)
    nA, nB, uA, uB, TA, TB = [Array{Float64}(2) for i = 1:6]
    sgn1 = [0, 1]
    sgn2 = [1, 0]
	ε = 1.e-10
	vel = Array{Float64}(nz)
	mm = Array{Int64}(nedges)

	if smoothing
		#smooth out interface between nz0 and nz0+nsmooth
		#----------------------------------------------------
		if smooth_vars[1]==1
			println("density smoothing function")
			println("--------------------------")
			#for each ion species, find position of max and min gradients, then smooth out density
			for j = 1:nspec
				maxgrad = 0.
				mingrad = 1.e20
				for k = 1:nz-1
					grad = hydro.rho[k,j] / hydro.rho[k+1,j]
					if grad > maxgrad
						maxgrad = grad
						kgrad[1] = k
						nA[1] = hydro.rho[k,j] / mi[j]
						nB[1] = hydro.rho[k+1,j] / mi[j]
					end
					if grad < mingrad
						mingrad = grad
						kgrad[2] = k
						nA[2] = hydro.rho[k,j] / mi[j]
						nB[2] = hydro.rho[k+1,j] / mi[j]
					end
				end

				for k = 1:nedges
					if smooth_type[k,j]=="max"
						mm[k] = 1
					elseif smooth_type[k,j]=="min"
						mm[k] = 2
					end
				end

				for jj = 1:nedges
					m = mm[jj]
					for k = kgrad[m]-sgn1[m]*nsmooth : kgrad[m]+sgn2[m]*nsmooth
							#now smooth out interface
							y1 = log10(nA[m])
							y2 = log10(nB[m])
							x1 = log10(r[kgrad[m] - sgn1[m]*nsmooth]) #note r is in micron
							x2 = log10(r[kgrad[m] + sgn2[m]*nsmooth]) #note r is in micron
							b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
							a  = ( y1 - b ) / x1

							hydro.rho[k,j] = 10^(a*log10(r[k]) + b) * mi[j]
					end
				end
			end
		end

		if smooth_vars[2]==1
			println("velocity smoothing function")
			println("--------------------------")
			#for each ion species, find position of max and min gradients, then smooth out density
            ok = [0, 0]
			for j = 1:nspec
				maxgrad = 0.
				mingrad = 1.e20
				vel = max(ε, abs(hydro.u[:,j]))
				for k = 1:nz-1
					grad = vel[k] / vel[k+1]
					if grad > maxgrad
						maxgrad = grad
						kgrad[1] = k
						uA[1] = vel[k]
						uB[1] = vel[k+1]
						if (k+tsmooth<nz) & (k-tsmooth>0)
							ok[1] = 1
						end
					end
					if grad < mingrad
						mingrad = grad
						kgrad[2] = k
						uA[2] = vel[k]
						uB[2] = vel[k+1]
						if (k+tsmooth < nz) & (k-tsmooth > 0)
							ok[2] = 1
						end
					end
				end

				for k = 1:nedges
					if smooth_type[k,j]=="max"
						mm[k]=1
					elseif smooth_type[k,j]=="min"
						mm[k]=2
					end
				end

				for jj = 1:nedges
					m = mm[jj]
					if ok[m].ne.1
                        continue #because it is a bad point...
                    end

					for k = kgrad[m]-sgn1[m]*vsmooth : kgrad[m]+sgn2[m]*vsmooth
							#now smooth out interface
							y1 = log10(uA[m])
							y2 = log10(uB[m])
							x1 = log10(r[kgrad[m] - sgn1[m]*vsmooth]) #note r is in micron
							x2 = log10(r[kgrad[m] + sgn2[m]*vsmooth]) #note r is in micron
							b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
							a  = ( y1 - b ) / x1

							hydro.u[k,j] = - (10^( a*log10(r[k]) + b ))
					end

				end
			end
		end

		if smooth_vars[3]==1
			println("temperature smoothing function")
			println("------------------------------")
			#for each ion species, find position of max and min gradients, then smooth out density
			ok = [0, 0]
			for j = 1:nspec
				maxgrad = 0.
				mingrad = 1.e20
				for k = 1:nz-1
					grad = hydro.T[k,j] / hydro.T[k+1,j]
					if grad > maxgrad
						maxgrad = grad
						kgrad[1] = k
						TA[1] = hydro.T[k,j]
						TB[1] = hydro.T[k+1,j]
						if (k+tsmooth < nz) & (k-tsmooth > 0)
							ok[1] = 1
						end
					end
					if grad < mingrad
						mingrad = grad
						kgrad[2] = k
						TA[2] = hydro.T[k,j]
						TB[2] = hydro.T[k+1,j]
						if (k+tsmooth < nz) & (k-tsmooth > 0)
							ok[2] = 1
						end
					end
				end

				for k = 1:nedges
					if smooth_type[k,j]=="max"
						mm[k] = 1
					elseif smooth_type[k,j]=="min"
						mm[k] = 2
					end
				end

				for jj = 1:nedges
					m = mm[jj]
					if ok[m]!=1
                        continue #because it is a bad point...
                    end
					for k = kgrad[m]-sgn1[m]*tsmooth:kgrad[m] + sgn2[m]*tsmooth
							#now smooth out interface
							y1 = log10(TA[m])
							y2 = log10(TB[m])
							x1 = log10(r[kgrad[m] - sgn1[m]*tsmooth]) #note r is in micron
							x2 = log10(r[kgrad[m] + sgn2[m]*tsmooth]) #note r is in micron
							b  = ( y1*x2 - y2*x1 ) / ( x2 - x1 )
							a  = ( y1 - b ) / x1

							hydro.T[k,j] = 10^(a*log10(r[k]) + b)
					end

				end
			end
		end

		#now update pressure after all this smoothing...
		for j = 1:nspec
			hydro.p[:,j] = hydro.rho[:,j] / mi[j] .* hydro.T[:,j]
		end
	end
    return hydro
end

#-----------------------------------
function init_fluxes(hydro)

	println()
	println(" WARNING - forcing quasi-neutrality and zero-current conditions")
	println()

	hydro.rho[:,nspec+1] = 0.
	hydro.u[:,nspec+1] = 0.

	#electrons
    for k = 1:nz
    	for j = 1:nspec
			hydro.rho[k,nspec+1] = hydro.rho[k,nspec+1] + me * Zi[j] * hydro.rho[k,j] / mi[j] #quasi-neutrality
			hydro.u[k,nspec+1] = hydro.u[k,nspec+1] + Zi[j] * hydro.rho[k,j] * hydro.u[k,j] / mi[j]   #zero-current condition
    	end
    	if !restart #take species 1 as reference
			hydro.p[k,nspec+1] = hydro.rho[k,nspec+1] / me * hydro.T[k,1] #electron pressure
			hydro.T[k,nspec+1] = hydro.T[k,1] #electron temperature
		else #we have read the electron temperature from file
			hydro.p[k,nspec+1] = hydro.rho[k,nspec+1] / me * hydro.T[k,nspec+1] #electron pressure
		end
    end

#   correct for electron velocity
    hydro.u[:,nspec+1] = hydro.u[:,nspec+1] ./ (hydro.rho[:,nspec+1] / me)
#	apply BC at piston side
    for j = 1:nspec+1
        hydro.u[1, j] = 0.
    end

#	finally, calculate U1D and F1D for all species
	for j = 1:nspec
		hydro.U1D[:, 3*(j-1)+1] = hydro.rho[:,j]
		hydro.U1D[:, 3*(j-1)+2] = hydro.rho[:,j] .* hydro.u[:,j]
		hydro.U1D[:, 3*(j-1)+3] = hydro.p[:,j] / (g-1) + 0.5 * hydro.rho[:,j] .* hydro.u[:,j].^2
		hydro.F1D[:, 3*(j-1)+1] = hydro.rho[:,j] .* hydro.u[:,j]
		hydro.F1D[:, 3*(j-1)+2] = hydro.rho[:,j] .* hydro.u[:,j].^2 + hydro.p[:,j]
		hydro.F1D[:, 3*(j-1)+3] = hydro.u[:,j] .* ( g / (g-1) * hydro.p[:,j] + 0.5 * hydro.rho[:,j] .* hydro.u[:,j].^2 )
	end
	#electrons
	hydro.U1D[:,neqi+1] = hydro.p[:,nspec+1] ./ (g-1) + 0.5 * hydro.rho[:,nspec+1] .* hydro.u[:,nspec+1].^2
	hydro.F1D[:,neqi+1] = hydro.u[:,nspec+1] .* ( g / (g-1) * hydro.p[:,nspec+1] + 0.5 * hydro.rho[:,nspec+1] .* hydro.u[:,nspec+1].^2 )

    return hydro
end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
@testset "Benchmark systems over rationals" begin

    systems = [
        ("noon 5", Groebner.change_ordering(Groebner.noonn(5), :degrevlex)),
        ("noon 6", Groebner.change_ordering(Groebner.noonn(6), :degrevlex)),
        ("noon 7", Groebner.change_ordering(Groebner.noonn(7), :degrevlex)),
        ("cyclic 6", Groebner.change_ordering(Groebner.cyclicn(6), :degrevlex)),
        ("cyclic 7", Groebner.change_ordering(Groebner.cyclicn(7), :degrevlex)),
        ("katsura 5", Groebner.change_ordering(Groebner.katsuran(5), :degrevlex)),
        ("katsura 6", Groebner.change_ordering(Groebner.katsuran(6), :degrevlex)),
        ("katsura 7", Groebner.change_ordering(Groebner.katsuran(7), :degrevlex)),
        ("henrion 5", Groebner.change_ordering(Groebner.henrion5(), :degrevlex)),
        ("henrion 6", Groebner.change_ordering(Groebner.henrion6(), :degrevlex)),
    ]

    for (name, system) in systems
        @info "testing" name
        @test Groebner.isgroebner(Groebner.groebner(system))
    end

end


text:
#=

    This is a port of the R2 model Grass.cs

    We observe that the grass is wet. Did it rain?

    Output from the R2 model:
    ```
    Mean: 0.699
    Variance: 0.21061
    Number of accepted samples = 790
    ```

    This model:

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

      cloudy    0.5530    0.4974     0.0157    0.0225   518.9849    0.9997       63.1061
        rain    0.6890    0.4631     0.0146    0.0215   439.4618    1.0015       53.4365
   sprinkler    0.4240    0.4944     0.0156    0.0196   422.0078    1.0014       51.3142
       temp1    0.7100    0.4540     0.0144    0.0197   523.3179    0.9990       63.6330
     wetRoof    0.4850    0.5000     0.0158    0.0231   457.4633    0.9992       55.6254
       temp2    0.9690    0.1734     0.0055    0.0107   514.7809    1.0017       62.5950
       temp3    0.9280    0.2586     0.0082    0.0125   565.5121    0.9990       68.7636
    wetGrass    1.0000    0.0000     0.0000    0.0000        NaN       NaN           NaN

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

# Did it rain?
@model function grass()
    cloudy ~ flip(0.5)

    if cloudy
        rain ~ flip(0.8)
        sprinkler ~ flip(0.1)
    else
        rain ~ flip(0.2)
        sprinkler ~ flip(0.5)
    end

    temp1 ~ flip(0.7)

    wetRoof ~ Dirac(temp1 && rain)

    temp2 ~ flip(0.9)
    temp3 ~ flip(0.9)

    wetGrass ~ Dirac(temp2 && rain || temp3 && sprinkler)

    # We observe that the grass is wet.
    # Did it rain?
    true ~ Dirac(wetGrass)

end

model = grass()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 1_000)

# chns = sample(model, NUTS(), 1_000)
# chns = sample(model, HMC(0.1,6), 1_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using OperatorLearning
using Documenter, DocumenterTools

DocMeta.setdocmeta!(OperatorLearning, :DocTestSetup, :(using OperatorLearning); recursive=true)

makedocs(;
    modules=[OperatorLearning],
    authors="Patrick Zimbrod <patrick.zimbrod@gmail.com> and contributors",
    repo="https://github.com/pzimbrod/OperatorLearning.jl/blob/{commit}{path}#{line}",
    sitename="OperatorLearning.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pzimbrod.github.io/OperatorLearning.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" =>
            ["Burgers Equation with FNO" => "examples/burgers_FNO.md",
             "Burgers Equation with DeepONet" => "examples/burgers_DeepONet.md"],
        "Frequently Asked Questions" => "faq.md",
        "Module Reference" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/pzimbrod/OperatorLearning.jl",
    #devbranch="main",
)


text:
#!/usr/bin/env julia

packages = [
    "Plots",
    "PyPlot",
    "PlotlyJS",
    "IJulia"
]

Pkg.update()

for p in packages
    @printf "... Install %s\n" p
    Pkg.add(p)
end

withenv("PYTHON" => "") do
    Pkg.build("PyPlot")
    Pkg.build("PlotlyJS")
end

try
    @printf "... Install AtomShell\n"
    using Blink
    Blink.AtomShell.install()
catch e
    @printf "Install failed: AtomShell\n"
    rethrow(e)
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
        "--min_nv"
            help = "Minimum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemin(Int)
        "--max_nv"
            help = "Maximum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemax(Int)
        "--recompute"
            help = "Wether to recompute the basic features if they're already not NULL in the database."
            default = false
            action = :store_true
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    db = SQLite.DB(parsed_args["dbpath"])
    save_features_instances!(db;
                             min_nv=parsed_args["min_nv"], max_nv=parsed_args["max_nv"],
                             recompute=parsed_args["recompute"])
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
    function ismv(
        f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)},
        mask::AbstractArray{Bool, 3},
        vsz::NTuple{3, Real};
        r::Real = 2*maximum(vsz),
        tol::Real = 1e-3,
        maxit::Integer = 500,
        verbose::Bool = false,
    ) -> Tuple{typeof(similar(f)), typeof(similar(mask))}

Iterative spherical mean value method (iSMV) [1].

### Arguments
- `f::AbstractArray{<:AbstractFloat, N ∈ (3, 4)}`: unwrapped (multi-echo) field/phase
- `mask::AbstractArray{Bool, 3}`: binary mask of region of interest
- `vsz::NTuple{3, Real}`: voxel size for smv kernel

### Keywords
- `r::Real = 2*maximum(vsz)`: radius of smv kernel in units of `vsz`
- `tol::Real = 1e-3`: stopping tolerance
- `maxit::Integer = 500`: maximum number of iterations
- `verbose::Bool = false`: print convergence information

### Returns
- `typeof(similar(f))`: background corrected local field/phase
- `typeof(similar(mask))`: eroded binary mask

### References
[1] Wen Y, Zhou D, Liu T, Spincemaille P, Wang Y. An iterative spherical mean
    value method for background field removal in MRI.  Magnetic resonance in
    medicine. 2014 Oct;72(4):1065-71.
"""
function ismv(
    f::AbstractArray{<:AbstractFloat, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real};
    r::Real = 2*maximum(vsz),
    tol::Real = 1e-3,
    maxit::Integer = 500,
    verbose::Bool = false,
) where {N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))
    return _ismv!(tzero(f), tzero(mask), f, mask, vsz, r, tol, maxit, verbose)
end

function _ismv!(
    fl::AbstractArray{<:AbstractFloat, N},
    smask::AbstractArray{Bool, 3},
    f::AbstractArray{T, N},
    mask::AbstractArray{Bool, 3},
    vsz::NTuple{3, Real},
    r::Real,
    tol::Real,
    maxit::Integer,
    verbose::Bool,
) where {T<:AbstractFloat, N}
    N ∈ (3, 4) || throw(ArgumentError("arrays must be 3d or 4d, got $(N)d"))

    size(fl) == size(f)        || throw(DimensionMismatch())
    size(smask) == size(mask)  || throw(DimensionMismatch())
    size(mask) == size(f)[1:3] || throw(DimensionMismatch())

    # crop image and pad for convolution
    Rc = crop_indices(mask)
    pad = ntuple(_ -> Int(cld(r, minimum(vsz))), Val(3))

    fp = padfastfft(@view(f[Rc,1]), pad, rfft=true)
    m = padfastfft(@view(mask[Rc]), pad, rfft=true)

    # init vars and fft plans
    sz = size(m)
    sz_ = (sz[1]>>1 + 1, sz[2], sz[3])

    s = similar(fp)
    fb = similar(fp)
    bc = similar(fp)
    m0 = similar(m)

    S = Array{T}(undef, sz_)
    F̂ = Array{complex(T)}(undef, sz_)

    FFTW.set_num_threads(num_cores())
    P = plan_rfft(fp)
    iP = inv(P)

    # get smv kernel
    S = _smv_kernel!(S, F̂, s, vsz, r, P)

    # constants
    _δ = one(eltype(s)) - sqrt(eps(eltype(s)))

    # erode mask
    s = _tcopyto!(s, m) # in-place type conversion, reuse smv var
    m0 = _tcopyto!(m0, m)

    F̂ = mul!(F̂, P, s)
    @inbounds @batch for I in eachindex(F̂)
        F̂[I] *= S[I]
    end

    s = mul!(s, iP, F̂)
    @inbounds @batch for I in eachindex(m)
        m[I] = s[I] > _δ
    end

    # iSMV
    @inbounds for t in axes(f, 4)
        if verbose && size(f, 4) > 1
            @printf("Echo: %d/%d\n", t, size(f, 4))
        end

        if t > 1
            fp = padarray!(fp, @view(f[Rc,t]))
        end

        @batch for I in eachindex(bc)
            bc[I] = (m0[I] - m[I]) * fp[I]
            fb[I] = fp[I]
        end

        fb = __ismv!(fb, s, S, bc, m, iP, F̂, P, tol, maxit, verbose)

        @batch for I in eachindex(fp)
            fp[I] = m[I] * (fp[I] - fb[I])
        end

        unpadarray!(@view(fl[Rc,t]), fp)
    end

    unpadarray!(@view(smask[Rc]), m)

    return fl, smask
end

function __ismv!(f::AbstractArray{T}, f0, S, bc, m, iP, F̂, P, tol, maxit, verbose) where {T}
    @inbounds @batch for I in eachindex(f0)
        f0[I] = m[I]*f[I]
    end

    nr = norm(f0)
    ϵ = tol * nr

    if verbose
        @printf("iter%6s\tresnorm\n", "")
    end

    @inbounds for i in 1:maxit
        if nr ≤ ϵ
            break
        end

        F̂ = mul!(F̂, P, f)
        @batch for I in eachindex(F̂)
            F̂[I] *= S[I]
        end

        f = mul!(f, iP, F̂)
        @batch threadlocal=zero(T)::T for I in eachindex(f)
            f[I] = muladd(m[I], f[I], bc[I])

            # norm(f0 - f)
            f0_f = f0[I] - f[I]
            threadlocal = muladd(f0_f, f0_f, threadlocal)
            f0[I] = f[I]
        end

        nr = sqrt(sum(threadlocal::Vector{T}))

        if verbose
            @printf("%5d/%d\t%1.3e\n", i, maxit, nr)
        end
    end

    if verbose
        println()
    end

    return f
end


text:
# type 1
# ------------------------------------------------------------------------------

const GAE{C,V} = JuMP.GenericAffExpr{C,V}
const GAEv{C} = JuMP.GenericAffExpr{C,JuMP.VariableRef}
const GAEp{C} = JuMP.GenericAffExpr{C,ParameterRef}

mutable struct DoubleGenericAffExpr{C,V,P} <: JuMP.AbstractJuMPScalar
    v::JuMP.GenericAffExpr{C,V}
    p::JuMP.GenericAffExpr{C,P}
end
const DGAE{C,V,P} = DoubleGenericAffExpr{C,V,P}

const ParametrizedGenericAffExpr{C,V} = DoubleGenericAffExpr{C, V, ParameterRef}
const PGAE{C,V} = ParametrizedGenericAffExpr{C,V}

const ParametrizedAffExpr{C} = DoubleGenericAffExpr{C, JuMP.VariableRef, ParameterRef}
const PAE{C} = ParametrizedAffExpr{C}
const PAEC{S} = JuMP.ScalarConstraint{PAE{Float64}, S}
const PVAEC{S} = JuMP.VectorConstraint{PAE{Float64}, S}

Base.one(::Type{ParameterRef}) = one(GAEp{Float64})

Base.iszero(a::PAE) = iszero(a.v) && iszero(a.p)
Base.zero(::Type{DGAE{C,V,P}}) where {C,V,P} = DGAE{C,V,P}(zero(GAEv{C}), zero(GAEp{C}))
Base.one(::Type{DGAE{C,V,P}}) where {C,V,P} = DGAE{C,V,P}(one(GAEv{C}), zero(GAEp{C}))
Base.zero(a::PAE) = zero(typeof(a))
Base.one( a::PAE) =  one(typeof(a))
Base.copy(a::DGAE{C,V,P}) where {C,V,P}  = DGAE{C,V,P}(copy(a.v), copy(a.p))
Base.broadcastable(expr::PAE) = Ref(expr)
# JuMP.GenericAffExpr{C,ParameterRef}(params::Vector{ParameterRef},coefs::Vector{C}) = JuMP.GenericAffExpr{C}(params,coefs,C(0.0))

DGAE{C,V,P}() where {C,V,P} = zero(DGAE{C,V,P})

Base.convert(::Type{PAE{C}}, aff::GAEv{C}) where {C} = PAE{C}(aff, GAEp{C}(zero(C)))

function Base.convert(::Type{PAE{T}}, aff::PAE) where {T}
    return PAE{T}(
        GAEv{T}(
            T(aff.v.constant),
            Pair{VariableRef,T}[k => T(v) for (k, v) in aff.v.terms],
        ),
        GAEp{T}(
            T(aff.p.constant),
            Pair{ParameterRef,T}[k => T(v) for (k, v) in aff.p.terms],
        ),
    )
end

function JuMP.map_coefficients_inplace!(f::Function, a::PAE)
    map_coefficients_inplace!(f, a.v)
    # The iterator remains valid if existing elements are updated.
    for (coef, par) in linear_terms(a.p)
        a.p.terms[par] = f(coef)
    end
    return a
end

function JuMP.map_coefficients(f::Function, a::PAE)
    return JuMP.map_coefficients_inplace!(f, copy(a))
end

function Base.sizehint!(a::PAE, n::Int, m::Int)
    sizehint!(a.v.terms, n)
    sizehint!(a.p.terms, m)
end

function JuMP.value(var_value::Function, ex::PAE{T})::T where {T}
    ret = value(var_value, ex.v)
    for (param, coef) in ex.p.terms
        ret += coef * var_value(param)
    end
    ret
end

JuMP.constant(aff::PAE) = aff.v.constant

# iterator

# add to expression

function Base.isequal(aff::DGAE{C,V,P}, other::DGAE{C,V,P}) where {C,V,P}
    return isequal(aff.v, other.v) && isequal(aff.p, other.p)
end

Base.hash(aff::DGAE, h::UInt) = hash(aff.v.constant, hash(aff.v.terms, h), hash(aff.p.terms, h))

function SparseArrays.dropzeros(aff::DGAE{C,V,P}) where {C,V,P}
    v = SparseArrays.dropzeros(aff.v)
    p = SparseArrays.dropzeros(aff.p)
    return DGAE{C,V,P}(v,p)
end

# Check if two AffExprs are equal after dropping zeros and disregarding the
# order. Mostly useful for testing.
function JuMP.isequal_canonical(aff::DGAE{C,V,P}, other::DGAE{C,V,P}) where {C,V,P}
    aff_nozeros = dropzeros(aff)
    other_nozeros = dropzeros(other)
    # Note: This depends on equality of OrderedDicts ignoring order.
    # This is the current behavior, but it seems questionable.
    return isequal(aff_nozeros, other_nozeros)
end
function Base.isapprox(x::GAEp, y::GAEp; kws...)
    if !isapprox(x.constant, y.constant; kws...)
        return false
    end
    x = dropzeros(x)
    y = dropzeros(y)
    if length(linear_terms(x)) != length(linear_terms(y))
        return false
    end
    for (coef, var) in linear_terms(x)
        c = get(y.terms, var, nothing)
        if c === nothing
            return false
        elseif !isapprox(coef, c; kws...)
            return false
        end
    end
    return true
end
function Base.isapprox(x::GAEp, y::PAE; kws...)
    return isapprox(x, y.p + y.v.constant; kws...)
end
function Base.isapprox(x::PAE, y::GAEp; kws...)
    return isapprox(y, x; kws...)
end

Base.convert(::Type{PAE{C}}, v::JuMP.VariableRef) where {C} = PAE{C}(GAEv{C}(zero(C), v => one(C)), zero(GAEp{C}))
Base.convert(::Type{PAE{C}}, p::ParameterRef) where {C} = PAE{C}(zero(GAEv{C}), GAEp{C}(zero(C), p => one(C)))
Base.convert(::Type{PAE{C}}, r::Real) where {C} = PAE{C}(GAEv{C}(convert(C, r)), zero(GAEp{C}))

# Check all coefficients are finite, i.e. not NaN, not Inf, not -Inf
function JuMP._assert_isfinite(a::DGAE)
    _assert_isfinite(v)
    for (coef, par) in linear_terms(a.p)
        isfinite(coef) || error("Invalid coefficient $coef on parameter $par.")
    end
end

JuMP.value(a::DGAE) = JuMP.value(value, a)

function JuMP.check_belongs_to_model(a::DGAE, model::AbstractModel)
    JuMP.check_belongs_to_model(a.v, model)
    JuMP.check_belongs_to_model(a.p, model)
end

# Note: No validation is performed that the variables in the AffExpr belong to
# the same model. The verification is done in `check_belongs_to_model` which
# should be called before calling `MOI.ScalarAffineFunction`.
# function MOI.ScalarAffineFunction(a::AffExpr)
#     _assert_isfinite(a)
#     terms = MOI.ScalarAffineTerm{Float64}[MOI.ScalarAffineTerm(t[1],
#                                                                index(t[2]))
#                                           for t in linear_terms(a)]
#     return MOI.ScalarAffineFunction(terms, a.constant)
# end
# moi_function(a::GenericAffExpr) = MOI.ScalarAffineFunction(a)
# function moi_function_type(::Type{<:GenericAffExpr{T}}) where T
#     return MOI.ScalarAffineFunction{T}
# end

# Copy an affine expression to a new model by converting all the
# variables to the new model's variables
# TODO function Base.copy(a::GenericAffExpr, new_model::AbstractModel)

# Build constraint
# ------------------------------------------------------------------------------

function JuMP.build_constraint(
    ::Function,
    aff::PAE{Float64},
    set::Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo},
)
    shifted_set = MOIU.shift_constant(set, -aff.v.constant)
    aff.v.constant = 0.0
    return JuMP.ScalarConstraint(aff, shifted_set)
end

function JuMP.build_constraint(
    errf::Function,
    aff::PAE,
    set::Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo},
)
    return build_constraint(errf, convert(PAE{Float64}, aff), set)
end

function JuMP.add_constraint(m::JuMP.Model, c::PAEC{S}, name::String="") where S
    # build LinearConstraint
    c_lin = JuMP.ScalarConstraint(c.func.v, c.set)
    # JuMP´s standard add_constrint
    cref = JuMP.add_constraint(m, c_lin, name)
    data = _getparamdata(m)::_ParameterData
    # needed for lazy get dual
    if lazy_duals(data)
        for (p, coef) in c.func.p.terms
            push!(data.constraints_map[index(p)], ParametrizedConstraintRef(cref, coef))
        end
    end
    # save the parameter part of a parametric affine expression
    _get_param_dict(data, S)[cref] = c.func.p
    return cref
end

# update constraint
# ------------------------------------------------------------------------------

function _update_constraints(data::_ParameterData)
    _update_constraints(data, MOI.EqualTo{Float64})
    _update_constraints(data, MOI.LessThan{Float64})
    _update_constraints(data, MOI.GreaterThan{Float64})
end

function _update_constraints(data::_ParameterData, ::Type{S}) where S
    for (cref, gaep) in _get_param_dict(data, S)
        _update_constraint(data, cref, gaep)
    end
    return nothing
end

function _update_constraint(data::_ParameterData, cref, gaep::GAEp{C}) where C
    val = 0.0
    @inbounds for (p, coef) in gaep.terms
        ind = index(p)
        val += coef * (data.future_values[ind] - data.current_values[ind])
    end
    _update_constraint(data, cref, val)
    return nothing
end

function _update_constraint(data::_ParameterData, cref, val::Number)
    if !iszero(val)
        JuMP.add_to_function_constant(cref, val)
    end
    return nothing
end

# Duals
# ------------------------------------------------------------------------------

function JuMP.dual(p::ParameterRef)
    params = _getparamdata(p)::_ParameterData
    if lazy_duals(params)
        return _getdual(p)
    else
        return params.dual_values[index(p)]
    end
end

# lazy

function _getdual(p::ParameterRef)
    return _getdual(p.model, index(p))
end
function _getdual(pcr::ParametrizedConstraintRef)::Float64
    pcr.coef * JuMP.dual(pcr.cref)
end
function _getdual(model::JuMP.Model, ind::Integer)
    params = _getparamdata(model)::_ParameterData
    # See (12) in http://www.juliaopt.org/MathOptInterface.jl/stable/apimanual.html#Duals-1
    # in the dual objective: - sum_i b_i^T y_i
    # Here b_i depends on the param and is b_i' + coef * param
    # so the dualobjective is:
    # - sum_i b_i'^T y_i - param * sum_i coef^T y_i
    # The dual of the parameter is: - sum_i coef^T y_i
    return -sum(_getdual.(params.constraints_map[ind]))
end

# non lazy

function _update_duals(data::_ParameterData)
    if !lazy_duals(data)
        fill!(data.dual_values, 0.0)
        _update_duals(data, MOI.EqualTo{Float64})
        _update_duals(data, MOI.GreaterThan{Float64})
        _update_duals(data, MOI.LessThan{Float64})
    end
    return nothing
end

function _update_duals(data::_ParameterData, ::Type{S}) where S
    for (cref, gaep) in _get_param_dict(data, S)
        _update_duals(data, cref, gaep)
    end
    return nothing
end

function _update_duals(data::_ParameterData, cref, gaep)
    dual_sol = JuMP.dual(cref)
    @inbounds for (p, coef) in gaep.terms
        data.dual_values[index(p)] -= coef * dual_sol
    end
    return nothing
end

# constraint modification
# ------------------------------------------------------------------------------

function JuMP.set_normalized_coefficient(con::CtrRef{F, S}, p::ParameterRef, coef::Number) where {F<:SAF, S}
    data = _getparamdata(p)
    dict = _get_param_dict(data, S)
    old_coef = 0.0
    ind = index(p)
    if haskey(dict, con)
        gaep = dict[con]
        old_coef = get!(gaep.terms, p, 0.0)
        gaep.terms[p] = coef
    else
        # TODO fix type C
        dict[con] = GAEp{Float64}(zero(Float64), p => coef)
    end
    if !iszero(coef-old_coef)
        val = (coef-old_coef)*data.current_values[ind]
        _update_constraint(data, con, val)
        if !iszero(data.future_values[ind] - data.current_values[ind])
            data.sync = false
        end
    end
    if lazy_duals(data)
        ctr_map = data.constraints_map[ind]
        found_ctr = false
        if isempty(ctr_map)
            push!(ctr_map, ParametrizedConstraintRef(con, coef))
        else
            for (p_ind, pctr) in enumerate(ctr_map)
                if pctr.cref == con
                    if found_ctr
                        deleteat!(ctr_map, p_ind)
                    else
                        ctr_map[p_ind] = ParametrizedConstraintRef(con, coef)
                    end
                    found_ctr = true
                end
            end
        end
        if found_ctr && !iszero(data.future_values[ind])
            data.sync = false
        end
    end
    nothing
end

"""
    delete_from_constraint(con, param::ParameterRef)

Removes parameter `param` from constraint `con`.
"""
function delete_from_constraint(con::CtrRef{F, S}, p::ParameterRef) where {F, S}
    data = _getparamdata(p)
    dict = _get_param_dict(data, S)
    ind = index(p)
    if haskey(dict, con)
        old_coef = get!(dict[con].terms, p, 0.0)
        _update_constraint(data, con, (0.0-old_coef) * data.current_values[ind])
        delete!(dict[con].terms, p)
        if !iszero(data.future_values[ind] - data.current_values[ind])
            data.sync = false
        end
    end
    if lazy_duals(data)
        ctr_map = data.constraints_map[ind]
        found_ctr = false
        for (index, pctr) in enumerate(ctr_map)
            if pctr.cref == con
                deleteat!(ctr_map, index)
                found_ctr = true
            end
        end
        if found_ctr && !iszero(data.future_values[ind])
            data.sync = false
        end
    end
    nothing
end

"""
    delete_from_constraints(param::ParameterRef)

Removes parameter `param` from all constraints.
"""
function delete_from_constraints(::Type{S}, p::ParameterRef) where S
    data = _getparamdata(p)
    dict = _get_param_dict(data, S)
    ind = index(p)
    for (con, gaep) in dict
        if haskey(gaep.terms, p)
            if !iszero(gaep.terms[p]) && !iszero(data.future_values[ind])
                data.sync = false
            end
            old_coef = get!(dict[con].terms, p, 0.0)
            _update_constraint(data, con, (0.0-old_coef) * data.current_values[ind])
            delete!(gaep.terms, p)
        end
    end
    if lazy_duals(data)
        if !isempty(data.constraints_map[ind])
            data.constraints_map[ind] = ParametrizedConstraintRef[]
            if !iszero(data.future_values[ind])
                data.sync = false
            end
        end
    end
    nothing
end

function delete_from_constraints(param::ParameterRef)
    delete_from_constraints(MOI.EqualTo{Float64}, param)
    delete_from_constraints(MOI.LessThan{Float64}, param)
    delete_from_constraints(MOI.GreaterThan{Float64}, param)
    # TODO lazy
    nothing
end

max_stars_repo_name: joaquimg/ParameterJuMP

max_issues_repo_name: joaquimg/ParameterJuMP

max_forks_repo_name: joaquimg/ParameterJuMP


text:
module BaseTraits
using SimpleTraits

export IsLeafType, IsBits, IsImmutable, IsContiguous, IsFastLinearIndex,
       IsAnything, IsNothing, IsCallable

"Trait which contains all types"
@traitdef IsAnything{X}
SimpleTraits.trait{X}(::Type{IsAnything{X}}) = IsAnything{X}

"Trait which contains no types"
typealias IsNothing{X} Not{IsAnything{X}}

"Trait of all isbits-types"
@traitdef IsBits{X}
@generated SimpleTraits.trait{X}(::Type{IsBits{X}}) =
    isbits(X) ? :(IsBits{X}) : :(Not{IsBits{X}})

"Trait of all immutable types"
@traitdef IsImmutable{X}
@generated SimpleTraits.trait{X}(::Type{IsImmutable{X}}) =
    X.mutable ? :(Not{IsImmutable{X}}) : :(IsImmutable{X})

"Trait of all callable objects"
@traitdef IsCallable{X}
@generated SimpleTraits.trait{X}(::Type{IsCallable{X}}) =
    (X==Function ||  length(methods(call, (X,Vararg)))>0) ? IsCallable{X} : Not{IsCallable{X}}

"Trait of all leaf types types"
@traitdef IsLeafType{X}
@generated trait{X}(::Type{IsLeafType{X}}) = isleaftype(X) ? :(Not{IsLeafType{X}}) : :(IsLeafType{X})

"Types which have contiguous memory layout"
@traitdef IsContiguous{X} # https://github.com/JuliaLang/julia/issues/10889
@generated SimpleTraits.trait{X}(::Type{IsContiguous{X}}) =
    Base.iscontiguous(X) ? :(IsContiguous{X}) : :(Not{IsContiguous{X}})

"Array indexing trait."
@traitdef IsFastLinearIndex{X} # https://github.com/JuliaLang/julia/pull/8432
@generated function SimpleTraits.trait{X}(::Type{IsFastLinearIndex{X}})
    if Base.linearindexing(X)==Base.LinearFast()
        return :(IsFastLinearIndex{X})
    elseif  Base.linearindexing(X)==Base.LinearSlow()
        return :(Not{IsFastLinearIndex{X}})
    else
        error("Not recognized")
    end
end

# TODO
## @traitdef IsArray{X} # use for any array like type in the sense of container
##                   # types<:AbstractArray are automatically part

## @traitdef IsMartix{X} # use for any LinearOperator
##                    # types<:AbstractArray are automatically part

end # module


text:
module NISTCrossSectionData
using ArgCheck
using Unitful
const UF = Unitful
using Unitful: MeV, g, cm

const uenergy  = MeV
const umassatt = cm^2/g
const umassrange = inv(umassatt)
const ustoppow = uenergy * umassatt
const unone    = Unitful.NoUnits

export MeV, g, cm

include("api.jl")
include("tables.jl")
include("csv.jl")
include("xaamdi.jl")
include("xcom.jl")
include("estar.jl")

end # module


text:
module orthoballoon_examples
using FinEtools

function orthoballoon()
    # Orthotropic balloon inflation, axially symmetric model

    # Parameters:
    E1=1.0;
    E2=1.0;
    E3=3.0;
    nu12=0.29;
    nu13=0.29;
    nu23=0.19;
    G12=0.3;
    G13=0.3;
    G23=0.3;
    p= 0.15;
    rin=1.;
    rex =1.2;
    tolerance=rin/1000.

    ##
    # Note that the FinEtools objects needs to be created with the proper
    # model-dimension reduction at hand.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DAxisymm

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the thickness in the radial direction, and the second
    # coordinate is the thickness in the axial direction.

    fens,fes = Q4block(rex-rin,pi/2,5,20);
    bdryfes = meshboundary(fes);

    icl = selectelem(fens, bdryfes, box=[0.,0.,0.,pi/2],inflate=tolerance);
    for i=1:count(fens)
        r=rin+fens.xyz[i,1]; a=fens.xyz[i,2];
        fens.xyz[i,:]=[r*cos(a) r*sin(a)];
    end

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # the symmetry plane
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    # the axis of symmetry
    l1 =selectnode(fens; box=[0 0 0 rex], inflate = tolerance)
    setebc!(u,l1,true, 1, 0.0)

    applyebc!(u)
    numberdofs!(u)
    println("Number of degrees of freedom = $(u.nfreedofs)")

    # The traction boundary condition is applied in the radial
    # direction.

    el1femm =  FEMMBase(IntegData(subset(bdryfes,icl), GaussRule(1, 3), true))
    function pressureloading!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copy!(forceout, XYZ/norm(XYZ)*p)
        return forceout
    end
    fi = ForceIntensity(FFlt, 2, pressureloading!); # pressure normal to the internal cylindrical surface
    F2= distribloads(el1femm, geom, u, fi, 2);

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2), true), material)

    K =stiffness(femm, geom, u)
    U=  K\(F2)
    scattersysvec!(u,U[:])

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 3)
    println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")

    File =  "orthoballoon_sigmaz.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
#pub_thick_pipe_axi()

end # orthoballoon

function orthoballoon_penalty()
    # Orthotropic balloon inflation, axially symmetric model

    # Parameters:
    E1=1.0;
    E2=1.0;
    E3=3.0;
    nu12=0.29;
    nu13=0.29;
    nu23=0.19;
    G12=0.3;
    G13=0.3;
    G23=0.3;
    p= 0.15;
    rin=1.;
    rex =1.2;
    tolerance=rin/1000.

    MR = DeforModelRed2DAxisymm

    fens,fes = Q4block(rex-rin,pi/2,5,20);
    bdryfes = meshboundary(fes);

    icl = selectelem(fens, bdryfes, box=[0.,0.,0.,pi/2], inflate=tolerance);
    for i=1:count(fens)
        r=rin+fens.xyz[i,1]; a=fens.xyz[i,2];
        fens.xyz[i,:]=[r*cos(a) r*sin(a)];
    end

    # now we create the geometry and displacement fields
    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # the symmetry plane
    ly =selectelem(fens, bdryfes; box=[0 rex 0 0], inflate = tolerance)
    # the axis of symmetry
    lx =selectelem(fens, bdryfes; box=[0 0 0 rex], inflate = tolerance)

    # No EBC
    applyebc!(u)
    numberdofs!(u)
    println("Number of degrees of freedom = $(u.nfreedofs)")

    # The traction boundary condition is applied in the radial
    # direction.

    el1femm =  FEMMBase(IntegData(subset(bdryfes,icl), GaussRule(1, 3), true))
    function pressureloading!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        copy!(forceout, XYZ/norm(XYZ)*p)
        return forceout
    end
    fi = ForceIntensity(FFlt, 2, pressureloading!); # pressure normal to the internal cylindrical surface
    F2= distribloads(el1femm, geom, u, fi, 2);

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2), true), material)

    ##
    # The restraints of the nodes on the bounding cross-sections in the direction
    # of the normal to the plane of the cross-section  in the
    # circumferential direction are introduced using a penalty formulation.
    # For that purpose we introduce  a finite element model machine for the
    # surface  finite elements on the cross-sections.
    springcoefficient =1.0e9/ (abs(p)/E1)
    xsfemm = FEMMDeforWinkler(IntegData(subset(bdryfes, lx), GaussRule(1, 3), true))
    ysfemm = FEMMDeforWinkler(IntegData(subset(bdryfes, ly), GaussRule(1, 3), true))
    H = surfacenormalspringstiffness(xsfemm,  geom, u, springcoefficient, SurfaceNormal(FFlt, 3)) +
    surfacenormalspringstiffness(ysfemm,  geom, u, springcoefficient, SurfaceNormal(FFlt, 3))
    K =stiffness(femm, geom, u)
    U=  (K + H)\(F2)
    scattersysvec!(u,U[:])

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 3)
    println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")

    File =  "orthoballoon_penalty_sigmaz.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    vectors=[("u", u.values)])
    @async run(`"paraview.exe" $File`)
#pub_thick_pipe_axi()

end # orthoballoon_penalty

function allrun()
    println("#####################################################")
    println("# orthoballoon ")
    orthoballoon()
    println("#####################################################")
    println("# orthoballoon_penalty ")
    orthoballoon_penalty()
    return true
end # function allrun

end # module orthoballoon_examples


text:
#! /usr/bin/julia

# Rosetta Code, Van der Corput sequence

function vdc{T<:Integer}(n::T, b::T)
    sum([d*float(b)^-i for (i, d) in enumerate(digits(n, b))])
end

for i in 2:9
    print("   Base ", i)
    for j in 0:9
        print(@sprintf(" %8.6f", vdc(j, i)))
    end
    println()
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Combinatorics

function part1()
	lines = readlines("inputD8.txt")
	easies = Int[2,4,3,7]
	counter = 0

	for line in lines
		afterpipe = split(line, "|")[2]
		words = split(afterpipe)
		counter += sum([length(word) in easies for word in words])
	end

	return counter
end

function part2()
	lines = readlines("inputD8.txt")
	letters_to_digit = Dict{String, Int}("abcefg" => 0, "abcdefg" => 8, "abdefg" => 6, "abcdfg" => 9, "cf" => 1, "abdfg" => 5, "bcdf" => 4, "acdeg" => 2, "acf" => 7, "acdfg" => 3)

	ordered = ["a", "b", "c", "d", "e", "f", "g"]
	perms = collect(permutations(ordered))
	total = 0
	for line in lines
		splitbypipe = split(line, "|")
		beforepipe = split.(split(splitbypipe[1]), "")
		for perm in perms
			mapping = Dict(ordered .=> perm)

			replaced = [[get(mapping, j, j) for j in b] for b in beforepipe]
			sorted = join.(sort.(replaced))

			if all(r in keys(letters_to_digit) for r in sorted)
				# wow, we did it
				afterpipe = split.(split(splitbypipe[2]), "")

				afterpipereplaced = [[get(mapping, j, j) for j in a] for a in afterpipe]
				afterpipesorted = join.(sort.(afterpipereplaced))

				final_number = parse(Int, join([letters_to_digit[i] for i in afterpipesorted]))
				total += final_number
				break
			end
		end
	end

	return total
end

println(part1())
println(part2())


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/IainNZ/Coverage.jl
#######################################################################

using Coverage, Base.Test

cd(Pkg.dir("Coverage")) do
    j = Coveralls.process_file(joinpath("test","data","Coverage.jl"))
    analyze_malloc(joinpath("test","data"))
end

srcname = joinpath("data","testparser.jl")
covname = srcname*".cov"
isfile(covname) && rm(covname)
cmdstr = "include(\"$srcname\"); using Base.Test; @test f2(2) == 4"
run(`julia --code-coverage=user -e $cmdstr`)
r = Coveralls.process_file(srcname)
# The next one is the correct one, but julia & JuliaParser don't insert a line number after the 1-line @doc -> test
# See https://github.com/JuliaLang/julia/issues/9663 (when this is fixed, can uncomment the next line on julia 0.4)
# target = [nothing, nothing, nothing, nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0]
target = [nothing, nothing, nothing, nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0]
@test r["coverage"][1:length(target)] == target

covtarget = (sum(x->x != nothing && x > 0, target), sum(x->x != nothing, target))
@test coverage_file(srcname) == covtarget
@test coverage_folder("data") != covtarget


text:
#
# 1D bound state energies for a QM particle in an infinite well and potential V(x)
#

## Definition of the orthonormal basis states
@everywhere function φ(i) return  function (x) sin(π*i*x) end end

## Coupling strength
@everywhere λ = 1.0*π^2/4.0

## Definition of the centered harmonic potential V(x)
@everywhere function V(x) return λ*(x-0.5)^2 end

## Multiply a potential to calculate <φ1|V|φ2>
@everywhere function Vφ(i) return  function (x) sin(π*i*x)*V(x) end end

## Definition of the inner product <φ1|φ2>
@everywhere function ∫(φ1::Function,φ2::Function)
    ψ(x)=φ1(x)*φ2(x)::Float64
    return 2*quadgk(ψ,0,1,abstol=1e-8)[1]::Float64
end

##  Infinite well energies
@everywhere function ɛ(i::Int,j::Int)
    return  i == j ? i^2 : 0
end

# The numerical Hamiltonian
function hamiltonian(n::Int)
  h = SharedArray(Float64,n^2)
  @sync for j=1:n
    @parallel for i=j:n
      h[i+(j-1)*n] = ɛ(i,j) + ∫(φ(i),Vφ(j))::Float64
     end
  end;
  H = reshape(h,n,n)
  return H+transpose(H-triu(H))
end

# The exact Hamiltonian
function hamiltonian_exact(n::Int)
  return Float64[ i==j? i^2 + (1/12)*λ*(1-6/(i*π)^2) : 0.25*(1+(-1)^(i+j))*(1/(i-j)^2-1/(i+j)^2) for i=1:n,j=1:n ]
end

## Checking orthonomality
function orthonormal(n::Int)
    v = SharedArray(Float64,n^2)
    @sync for j=1:n
        @parallel for i=j:n
            v[i+(j-1)*n] = ∫(φ(i),φ(j))::Float64
        end
    end;
    m = reshape(v,n,n)
    return m+transpose(m-triu(m))
end

function main(n::Int)
   @time println("norm of difference = ", norm(eigvals(hamiltonian(n))-eigvals(hamiltonian_exact(n))))
end

main(400)

max_stars_repo_name: lmelwyn/julia-demo

max_issues_repo_name: lmelwyn/julia-demo

max_forks_repo_name: lmelwyn/julia-demo


text:
"""
    rand_tangent([rng::AbstractRNG,] x)

Returns a randomly generated tangent vector appropriate for the primal value `x`.
"""
rand_tangent(x) = rand_tangent(Random.GLOBAL_RNG, x)

rand_tangent(rng::AbstractRNG, x::Symbol) = DoesNotExist()
rand_tangent(rng::AbstractRNG, x::AbstractChar) = DoesNotExist()
rand_tangent(rng::AbstractRNG, x::AbstractString) = DoesNotExist()

rand_tangent(rng::AbstractRNG, x::Integer) = DoesNotExist()

rand_tangent(rng::AbstractRNG, x::T) where {T<:Number} = randn(rng, T)

rand_tangent(rng::AbstractRNG, x::StridedArray) = rand_tangent.(Ref(rng), x)

function rand_tangent(rng::AbstractRNG, x::T) where {T<:Tuple}
    return Composite{T}(rand_tangent.(Ref(rng), x)...)
end

function rand_tangent(rng::AbstractRNG, xs::T) where {T<:NamedTuple}
    return Composite{T}(; map(x -> rand_tangent(rng, x), xs)...)
end

function rand_tangent(rng::AbstractRNG, x::T) where {T}
    if !isstructtype(T)
        throw(ArgumentError("Non-struct types are not supported by this fallback."))
    end

    field_names = fieldnames(T)
    if length(field_names) > 0
        tangents = map(field_names) do field_name
            rand_tangent(rng, getfield(x, field_name))
        end
        return Composite{T}(; NamedTuple{field_names}(tangents)...)
    else
        return NO_FIELDS
    end
end


text:
struct Heisenberg{D} <: AbstractHeisenberg{D}
    size::NTuple{D, Int}
    periodic::Bool
    Heisenberg(size::Int...; periodic::Bool) = new{length(size)}(size, periodic)
end

Base.size(model::Heisenberg) = model.size

heisenberg_ij(nbit::Int, i::Int, j::Int=i+1) = put(nbit, i=>X)*put(nbit, j=>X) + put(nbit, i=>Y)*put(nbit, j=>Y) + put(nbit, i=>Z)*put(nbit, j=>Z)
const heisenberg_term = repeat(2, X, 1:2) + repeat(2, Y, 1:2) + repeat(2, Z, 1:2)

function hamiltonian(model::AbstractHeisenberg)
    nbit = nspins(model)
    sum([x[3]*heisenberg_ij(nbit, x[1], x[2]) for x in get_bonds(model)])*0.25
end

function get_bonds(model::Heisenberg{2})
    m, n = model.size
    cis = LinearIndices(model.size)
    bonds = Tuple{Int, Int, Float64}[]
    for i=1:m, j=1:n
        (i!=m || model.periodic) && push!(bonds, (cis[i,j], cis[i%m+1,j], 1.0))
        (j!=n || model.periodic) && push!(bonds, (cis[i,j], cis[i,j%n+1], 1.0))
    end
    bonds
end

function get_bonds(model::AbstractModel{1})
    nbit, = model.size
    [(i, i%nbit+1, 1.0) for i in 1:(model.periodic ? nbit : nbit-1)]
end

"""
    energy(config, model::AbstractHeisenberg; nbatch) -> Float64

Ground state energy by sampling Quantum circuit.
The hamiltonian is limited to Heisenberg and J1J2 Type.
"""
function energy(qpeps::QPEPSMachine, model::AbstractHeisenberg)
    local eng = 0.0
    for basis in [X, Y, Z]
        mres = gensample(qpeps, basis)
        for (i,j,w) in get_bonds(model)
            eng += w*(1-2*mean(mres[i] .⊻ mres[j]))
        end
    end
    eng/=4
end


text:
function BlackScholes(CallPutFlag, S, X, T, r, v)

  d1 = (log(S / X) + (r + v^2 / 2) * T) / (v * sqrt(T))
  d2 = d1 - v * sqrt(T)

  if CallPutFlag == "c"
    BS = S * CND(d1) - X * exp(-r * T) * CND(d2)
  elseif CallPutFlag == "p"
    BS = X * exp(-r * T) * CND(-d2) - S * CND(-d1)
  end

end


text:
using Test

include("grains.jl")

@testset "beginning squares" begin
    @test on_square(1) == 1
    @test on_square(2) == 2
    @test on_square(3) == 4
    @test on_square(4) == 8
    @test on_square(16) == 32768
    @test on_square(32) == 2147483648
    @test total_after(1) == 1
    @test total_after(3) == on_square(1) + on_square(2) + on_square(3)
end

@testset "ending squares" begin
    @test total_after(32) < total_after(64)
    @test on_square(64) == 9223372036854775808
    @test total_after(64) == 18446744073709551615
end

@testset "Invalid values" begin
    @testset "Zero" begin
        @test_throws DomainError on_square(0)
        @test_throws DomainError total_after(0)
    end

    @testset "Negative" begin
        @test_throws DomainError on_square(-1)
        @test_throws DomainError total_after(-1)
    end

    @testset "Greater than 64" begin
        @test_throws DomainError on_square(65)
        @test_throws DomainError total_after(65)
    end
end


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

struct DefaultApi <: SwaggerApi
    client::Swagger.Client
end

"""

Checks whether a domain name in the cloudapp.azure.com zone is available for use.
Param: location::String (required)
Param: domainNameLabel::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DnsNameAvailabilityResult
"""
function _swaggerinternal_checkDnsNameAvailability(_api::DefaultApi, location::String, domainNameLabel::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DnsNameAvailabilityResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/locations/{location}/CheckDnsNameAvailability", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "domainNameLabel", domainNameLabel)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function checkDnsNameAvailability(_api::DefaultApi, location::String, domainNameLabel::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_checkDnsNameAvailability(_api, location, domainNameLabel, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function checkDnsNameAvailability(_api::DefaultApi, response_stream::Channel, location::String, domainNameLabel::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_checkDnsNameAvailability(_api, location, domainNameLabel, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export checkDnsNameAvailability


text:
__precompile__()
module StrainRecon
"""
D = ddx(n), 1D derivative operator
"""
function ddx(n)
	return spdiagm((fill(-1.0,n),fill(1.0,n)),(0,1),n,n+1)
end
export sdiag

"""
S = sdiag(s) builds sparse diagonal matrix
"""
function sdiag(a::Vector)
	n = length(a)
	i = collect(1:n+1) # colptr
	j = collect(1:n)   # rowval
	return SparseMatrixCSC(n,n,i,j,a)
end

# Here go all the packages required
using Combinatorics

include("StrainReconParam.jl")
include("findM.jl")
include("findw.jl")
include("misfit.jl")
include("iplp.jl")
include("bcd.jl")
include("solveMAP.jl")

include("randomw.jl")
include("integratePosterior.jl")
include("asqp.jl")
include("entropyM.jl")
include("srDistance.jl")
include("randomM.jl")
include("intmat2binmat.jl")

hasGurobi=false
try
	using Gurobi
	using JuMP
	hasGurobi = true
catch
end
if hasGurobi
	include("solveConvexMAP.jl")
end

export warmupStrainRecon
function warmupStrainRecon()
	(m,n) = (8,2)
	M     = rand(Bool,m,n)
	R     = sdiag(ones(m))
	atrue = rand(n); atrue/=sum(atrue)
	atrue = sort(atrue,rev=true)
	Data  = M*atrue

	param = getStrainReconParam(Data,n,R=R,label="init_test")

	println("getReduced Misfit")
	@time reducedMisfit(param,atrue; doDerivative=true)

	println("solveMAP")
    @time solveMAP(param;maxTry=1,storeInterm=false,maxIterQP=1,maxIterBCD=1)

end
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    abstract type SetMapBridge{T,S2,S1,F,G} <: AbstractBridge end

Consider two type of sets, `S1` and `S2`, and a linear mapping `A` such that
the image of a set of type `S1` under `A` is a set of type `S2`.

A `SetMapBridge{T,S2,S1,F,G}` is a bridge that maps `G`-in-`S2` constraints
into `F`-in-`S1` by mapping the function through `A`.

The linear map `A` is described by;

 * [`MathOptInterface.Bridges.map_set`](@ref)
 * [`MathOptInterface.Bridges.map_function`](@ref).

Implementing a method for these two functions is sufficient to bridge
constraints. However, in order for the getters and setters of attributes such as
dual solutions and starting values to work as well, a method for the following
functions must be implemented:

 * [`MathOptInterface.Bridges.inverse_map_set`](@ref)
 * [`MathOptInterface.Bridges.inverse_map_function`](@ref)
 * [`MathOptInterface.Bridges.adjoint_map_function`](@ref)
 * [`MathOptInterface.Bridges.inverse_adjoint_map_function`](@ref)

See the docstrings of each function to see which feature would be missing if it
was not implemented for a given bridge.
"""
abstract type SetMapBridge{T,S2,S1,F,G} <: AbstractBridge end

function bridge_constraint(
    BT::Type{<:SetMapBridge{T,S2,S1,F,G}},
    model::MOI.ModelLike,
    func::G,
    set::S1,
) where {T,S2,S1,F,G}
    mapped_func = MOI.Bridges.map_function(BT, func)
    mapped_set = MOI.Bridges.map_set(BT, set)
    constraint = MOI.add_constraint(model, mapped_func, mapped_set)
    return BT(constraint)
end

function MOI.supports_constraint(
    ::Type{<:SetMapBridge{T,S2,S1}},
    ::Type{<:MOI.AbstractScalarFunction},
    ::Type{S1},
) where {T,S2,S1<:MOI.AbstractScalarSet}
    return true
end

function MOI.supports_constraint(
    ::Type{<:SetMapBridge{T,S2,S1}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{S1},
) where {T,S2,S1<:MOI.AbstractVectorSet}
    return true
end

function MOI.Bridges.added_constrained_variable_types(::Type{<:SetMapBridge})
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:SetMapBridge{T,S2,S1,F}},
) where {T,S2,S1,F}
    return Tuple{Type,Type}[(F, S2)]
end

# Attributes, Bridge acting as a model

function MOI.get(
    ::SetMapBridge{T,S2,S1,F},
    ::MOI.NumberOfConstraints{F,S2},
)::Int64 where {T,S2,S1,F}
    return 1
end

function MOI.get(
    bridge::SetMapBridge{T,S2,S1,F},
    ::MOI.ListOfConstraintIndices{F,S2},
) where {T,S2,S1,F}
    return [bridge.constraint]
end

# References

function MOI.delete(model::MOI.ModelLike, bridge::SetMapBridge)
    MOI.delete(model, bridge.constraint)
    return
end

# Attributes, Bridge acting as a constraint

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::SetMapBridge{T,S2,S1,F,G},
) where {T,S2,S1,F,G}
    mapped_func = MOI.get(model, attr, bridge.constraint)
    func = MOI.Bridges.inverse_map_function(typeof(bridge), mapped_func)
    return MOI.Utilities.convert_approx(G, func)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::SetMapBridge{T,S2,S1,F,G},
    func::G,
) where {T,S2,S1,F,G}
    new_f = MOI.Bridges.map_function(typeof(bridge), func)
    MOI.set(model, attr, bridge.constraint, new_f)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SetMapBridge,
)
    set = MOI.get(model, attr, bridge.constraint)
    return MOI.Bridges.inverse_map_set(typeof(bridge), set)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SetMapBridge{T,S2,S1},
    set::S1,
) where {T,S2,S1}
    new_set = MOI.Bridges.map_set(typeof(bridge), set)
    MOI.set(model, attr, bridge.constraint, new_set)
    return
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:SetMapBridge{T,S2,S1,F}},
) where {T,S2,S1,F}
    return MOI.supports(model, attr, MOI.ConstraintIndex{F,S2})
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::SetMapBridge,
)
    value = MOI.get(model, attr, bridge.constraint)
    return MOI.Bridges.inverse_map_function(typeof(bridge), value)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::SetMapBridge,
    value,
)
    mapped_value = MOI.Bridges.map_function(typeof(bridge), value)
    MOI.set(model, attr, bridge.constraint, mapped_value)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::SetMapBridge,
)
    value = MOI.get(model, attr, bridge.constraint)
    return MOI.Bridges.adjoint_map_function(typeof(bridge), value)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::BT,
    value,
) where {BT<:SetMapBridge}
    mapped_value = MOI.Bridges.inverse_adjoint_map_function(BT, value)
    MOI.set(model, attr, bridge.constraint, mapped_value)
    return
end

# By linearity of the map, we can just change the constant/coefficient

function MOI.modify(
    model::MOI.ModelLike,
    bridge::BT,
    change::MOI.ScalarConstantChange,
) where {BT<:SetMapBridge}
    new_constant = MOI.Bridges.map_function(BT, change.new_constant)
    MOI.modify(model, bridge.constraint, MOI.ScalarConstantChange(new_constant))
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::BT,
    change::MOI.VectorConstantChange,
) where {BT<:SetMapBridge}
    new_constant = MOI.Bridges.map_function(BT, change.new_constant)
    MOI.modify(model, bridge.constraint, MOI.VectorConstantChange(new_constant))
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::BT,
    change::MOI.ScalarCoefficientChange,
) where {BT<:SetMapBridge}
    new_coefficient = MOI.Bridges.map_function(BT, change.new_coefficient)
    new_change = MOI.ScalarCoefficientChange(change.variable, new_coefficient)
    MOI.modify(model, bridge.constraint, new_change)
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::BT,
    change::MOI.MultirowChange,
) where {BT<:SetMapBridge}
    # It is important here that `change.new_coefficients` contains the complete
    # new sparse column associated to the variable. Calling modify twice with
    # part of the column won't work since  the linear map might reset all the
    # column each time.
    coefficients = MOI.Bridges.map_function(BT, change.new_coefficients)
    new_change = MOI.MultirowChange(change.variable, coefficients)
    MOI.modify(model, bridge.constraint, new_change)
    return
end


text:
module NLSolvers

import Base: show, summary

using Logging

#============================ LinearAlgebra ===========================
  We use often use the LinearAlgebra functions dot and norm for opera-
  tions related to assessing angles between vectors, size of vectors
  and so on.

  mul!, rmul!, ldiv!, ect can to make Array{T, N} operations as fast
  a possible.

  cholesky is impossible to live without to replace inverses

  UniformScaling (I), Symmetric, Diagonal are all usefull to handle
  shifted systems, hessians, etc

  eigen and diag are useful in trust region subprblems

  opnorm is used in the trust region subproblem safe guards to bound
  the estimate on lambda
 ============================ LinearAlgebra ===========================#

using LinearAlgebra:  dot, I, norm, # used everywhere in updates, convergence, etc
                      mul!, rmul!, ldiv!, # quasi-newton updates, apply factorizations, etc
                      cholesky, factorize, issuccess, # very useful in trust region solvers
                      UniformScaling, Diagonal, # simple matrices
                      Symmetric, Hermitian, # wrap before factorizations or eigensystems to avoid checks
                      diag, # mostly for trust region diagonal manipulation
                      eigen, # for the direct subproblem solver
                      opnorm, # for NWI safe guards
                      checksquare, UpperTriangular, givens, lmul!, cond, # For QR update
                      axpy! # for Anderson
import LinearAlgebra: mul!, dot # need to extend for preconditioners

# For better random number generators and rand!
using RandomNumbers

using Printf

function solve end
export solve
"""

"""
function objective_return(f, g, H=nothing)
  if g isa Nothing && H isa Nothing
    return f
  elseif !(g isa Nothing) && H isa Nothing
    return f, g
  elseif !(g isa Nothing) && !(H isa Nothing)
    return f, g, H
  end
end
export objective_return
isallfinite(x) = mapreduce(isfinite, *, x)

using StaticArrays
abstract type MutateStyle end

abstract type AbstractProblem end
abstract type AbstractOptions end

struct InPlace <: MutateStyle end
struct OutOfPlace <:MutateStyle end
include("precondition.jl")
include("Manifolds.jl")
include("objectives.jl")
include("linearalgebra.jl")
export NonDiffed, OnceDiffed, TwiceDiffed, ScalarObjective

# make this struct that has scheme and approx
abstract type QuasiNewton{T1} end

abstract type HessianApproximation end
struct Inverse <: HessianApproximation end
struct Direct <: HessianApproximation end
export Inverse, Direct
# problem and options types
include("optimize/problem_types.jl")
export OptimizationProblem, OptimizationOptions
# Globalization strategies
# TODO:

# Initial step
abstract type LineSearcher end
include("globalization/linesearches/root.jl")
export Backtracking, Static, HZAW
# step interpolations
export FFQuadInterp

include("globalization/trs_solvers/root.jl")
export NWI, Dogleg, NTR

# Quasi-Newton (including Newton and gradient descent) functionality
include("quasinewton/quasinewton.jl")
export DBFGS, BFGS, SR1, DFP, GradientDescent, Newton, BB, LBFGS, ActiveBox

# Include the actual functions that expose the functionality in this package.
include("optimize/linesearch/linesearch.jl")
export HS, CD, HZ, FR, PRP, VPRP, LS, DY
include("optimize/randomsearch/randomsearch.jl")
export SimulatedAnnealing, PureRandomSearch, APSO
include("optimize/directsearch/directsearch.jl")
export NelderMead
include("optimize/acceleration/root.jl")
export Adam, AdaMax

include("optimize/trustregions/trustregions.jl")
export minimize, minimize!, OptProblem, LineSearch, TrustRegion

include("optimize/projected/root.jl")

include("fixedpoints/root.jl")
export Anderson

include("nlsolve/root.jl")
export NEqProblem, NEqOptions
export InexactNewton, KrylovNEqProblem
export DFSANE
include("nlsolve/acceleration/anderson.jl")

# Forcing Terms
export FixedForceTerm, DemboSteihaug, EisenstatWalkerA, EisenstatWalkerB

function negate(problem::AbstractProblem, A)

end
function negate(::InPlace, A::AbstractArray)
  A .= .-A
  return A
end
function negate(::OutOfPlace, A::AbstractArray)
  -A
end

mstyle(problem::AbstractProblem) = problem.mstyle

"""
    retract(problem, z, x, p [, α])

Move from `x` along the direction `p` (or `α*p` if α is supplied) to a new point on the manifold in `problem` and store it in z if the problem is specified as inplace.  If the problem is inplace and updated `z` is returned, else a new vector is returned.
"""
retract(problem, z, x, p) = _retract(mstyle(problem), _manifold(problem), z, x, p)
retract(problem, z, x, p, α) = _retract(mstyle(problem), _manifold(problem), z, x, p, α)
function _retract(::InPlace, manifold::Manifold, z, x, p)
    retract!(manifold, z, x, p)
    return z
end
function _retract(::OutOfPlace, manifold::Manifold, z, x, p)
    z = retract(manifold, x, p)
    return z
end
# shouldn't have those here
function _retract(::InPlace, manifold::Euclidean, z, x, p, α)
    @. z = x + α*p
    return z
  end
function _retract(::OutOfPlace, manifold::Euclidean, z, x, p, α)
    z = @. x + α*p
    return z
end
function _retract(::InPlace, manifold::Euclidean, z, x, p)
    @. z = x + p
    return z
  end
function _retract(::OutOfPlace, manifold::Euclidean, z, x, p)
    z = @. x + p
    return z
end

"""
    project_tanget(problem, w, x, y)

Project g on the tangent space to the manifold (stored in `problem`) at x and store the result in `w` if the problem is specified as inplace. If the problem is inplace and updated `w` is returned, else a new vector is returned.
"""
function project_tangent(problem::OptimizationProblem, w, x, v)
    project_tanget(mstyle(problem), problem, w, x, v)
end
project_tanget(::InPlace, problem, w, x, v) = copyto!(w, v)
project_tanget(::OutOfPlace, problem, w, x, v) = v

end # module


text:
using Latexify
using Test

module MyModule
using Latexify
struct MyType
    vec1
    vec2
end

my_reverse(x) = x[end:-1:1]

@latexrecipe function f(t::MyType; reverse=false)
    env --> :align
    fmt := "%.2f"
    return reverse ? (t.vec2, t.vec1) : (t.vec1, t.vec2)
end

struct MyVec
    vec::AbstractVector
end

@latexrecipe function f(v::MyVec; reverse=false)
    env --> :array
    fmt := "%.2f"
    return reverse ? v.vec[end, -1, 1] : v.vec
end

struct MyTup
    tup::Tuple
end

@latexrecipe function f(v::MyTup; reverse=false)
    env --> :array
    fmt := "%.2f"
    return reverse ? my_reverse(v.tup) : v.tup
end

struct MyDoubleTup
    tup1::Tuple
    tup2::Tuple
end

@latexrecipe function f(t::MyDoubleTup)
    return t.tup1, t.tup2
end

struct MyFloat
    x::Float64
end

@latexrecipe function f(vec::Vector{T}) where T <: MyFloat
    fmt --> "%.4e"
    return [myfloat.x for myfloat in vec]
end

struct MyDoubleVec{T}
    vec1::AbstractVector{T}
    vec2::AbstractVector{T}
end

@latexrecipe function f(vec::MyDoubleVec{T}) where T <: MyFloat
    return [myfloat.x for myfloat in vec.vec1], [myfloat.x for myfloat in vec.vec2]
end

end

using .MyModule
t = MyModule.MyType([:A, :B, 3.], [1., 2, 3])
t2 = MyModule.MyType([:X, :Y, :(x/y)], Number[1.23434534, 232423.42345, 12//33])

vec = [MyModule.MyFloat(x) for x in 1:4]
@test latexify(vec; transpose=true) ==
raw"\begin{equation}
\left[
\begin{array}{cccc}
1.0000e+00 & 2.0000e+00 & 3.0000e+00 & 4.0000e+00 \\
\end{array}
\right]
\end{equation}
"

double_vec = MyModule.MyDoubleVec([MyModule.MyFloat(x) for x in 1:4], [MyModule.MyFloat(x) for x in 8:11])
@test latexify(double_vec) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1.0 & 8.0 \\
2.0 & 9.0 \\
3.0 & 10.0 \\
4.0 & 11.0 \\
\end{array}
\right]
\end{equation}
"

@test latexify(t2, fmt="%.8f") ==
raw"\begin{align}
X =& 1.23 \\
Y =& 232423.42 \\
\frac{x}{y} =& \frac{4.00}{11.00}
\end{align}
"

@test latexify(t) ==
raw"\begin{align}
A =& 1.00 \\
B =& 2.00 \\
3.00 =& 3.00
\end{align}
"

@test latexify(t; reverse=true) ==
raw"\begin{align}
1.00 =& A \\
2.00 =& B \\
3.00 =& 3.00
\end{align}
"

@test latexify(t; env=:array) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
A & 1.00 \\
B & 2.00 \\
3.00 & 3.00 \\
\end{array}
\right]
\end{equation}
"

@test latexify(t; env=:array, reverse=true) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1.00 & A \\
2.00 & B \\
3.00 & 3.00 \\
\end{array}
\right]
\end{equation}
"

vec = MyModule.MyVec([1., 2.])
@test latexify(vec, transpose=true) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1.00 & 2.00 \\
\end{array}
\right]
\end{equation}
"

tup = MyModule.MyTup((1., 2.))

@test latexify(tup, transpose=true) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1.00 & 2.00 \\
\end{array}
\right]
\end{equation}
"

@test latexify(tup, reverse=true, transpose=true) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
2.00 & 1.00 \\
\end{array}
\right]
\end{equation}
"

tup2 = MyModule.MyDoubleTup((1., 3), (2., 4.))
@test latexify(tup2) ==
raw"\begin{equation}
\left[
\begin{array}{cc}
1.0 & 2.0 \\
3 & 4.0 \\
\end{array}
\right]
\end{equation}
"


text:
using Test
using DFTK: load_psp, pymatgen_lattice, pymatgen_structure, load_lattice, load_atoms
using DFTK: ElementCoulomb
using Unitful
using UnitfulAtomic
using PyCall

py"""
from pymatgen.core.lattice import Lattice
from pymatgen.core.structure import Structure
"""

@testset "Lattice DFTK -> pymatgen -> DFTK" begin
    reference = randn(3, 3)
    output = load_lattice(pymatgen_lattice(reference))
    @test output ≈ reference atol=1e-14
end

@testset "Lattice pymatgen -> DFTK -> pymatgen" begin
    data = randn(9)
    reference = py"Lattice($data)"
    output = pymatgen_lattice(load_lattice(reference))
    output = py"$output.matrix.ravel()" .+ 0
    @test data ≈ output atol=1e-14
end

@testset "Lattice DFTK -> pymatgen" begin
    a = randn(3)
    b = randn(3)
    c = randn(3)
    reference = [a b c]

    # Convert the lattice to python, make it flat and convert
    # to Julia array
    output = pymatgen_lattice(reference)
    outlatt = py"$output.matrix.ravel()" .+ 0
    @test a ≈ austrip.(outlatt[1:3] * u"Å") atol=1e-14
    @test b ≈ austrip.(outlatt[4:6] * u"Å") atol=1e-14
    @test c ≈ austrip.(outlatt[7:9] * u"Å") atol=1e-14
end

@testset "Structure pymatgen -> DFTK -> pymatgen" begin
    data = randn(9)
    reflattice = py"Lattice($data)"
    species = [1, 1, 6, 6, 6, 8]
    positions = [randn(3), randn(3), randn(3), randn(3), randn(3), randn(3)]

    reference = py"Structure($reflattice, $species, $positions)"
    atoms = load_atoms(reference)
    @test length(atoms) == 3
    @test all(at isa ElementCoulomb for (at, positions) in atoms)
    @test atoms[1][1].symbol == :H
    @test atoms[2][1].symbol == :C
    @test atoms[3][1].symbol == :O

    output = pymatgen_structure(load_lattice(reference), atoms)
    @test output.lattice == reflattice
    for i in 1:6
        @test output.species[i].number == species[i]
        @test output.sites[i].frac_coords == positions[i]
    end
end

@testset "Structure DFTK -> pymatgen" begin
    a = randn(3)
    b = randn(3)
    c = randn(3)
    lattice = [a b c]

    atoms = [
        ElementCoulomb(1) => [randn(3), randn(3)],
        ElementCoulomb(6) => [randn(3), randn(3), randn(3)],
        ElementCoulomb(8) => [randn(3), randn(3)],
    ]

    # Convert the lattice to python, make it flat and convert
    # to Julia array
    output = pymatgen_structure(lattice, atoms)
    outlatt = py"$output.lattice.matrix.ravel()" .+ 0
    @test a ≈ austrip.(outlatt[1:3] * u"Å") atol=1e-14
    @test b ≈ austrip.(outlatt[4:6] * u"Å") atol=1e-14
    @test c ≈ austrip.(outlatt[7:9] * u"Å") atol=1e-14

    specmap = [1, 1, 2, 2, 2, 3, 3]
    offset = [0, 0, 2, 2, 2, 5, 5]
    for i in 1:6
        specpair = atoms[specmap[i]]
        ired = i - offset[i]
        @test output.species[i].number == specpair.first.Z
        @test output.sites[i].frac_coords == specpair.second[ired]
    end
end