text:
module GivensBoson
export given_eigen_solver
using LinearAlgebra

"""
    find the maximum-abs elements which is off-diagonal of ih.
    [If abs(a)<tol. a can be seen as zero.]
    return the value and the row and column index.
    and a. and If there is any off-diagonal element.

    Input manual:
        if zeromode = true: this function will save zero mode and avoid to return imformation about them.
        if zeromode = false: we suppose there is no zero mode in this problem.

    Output Format:value, index,flag

        value: ih[index...]
        flag:
            flag = 1: there is a non-zero-mode element we need to report.
            flag = 2: there is a zero mode we need to report.
            flag = 3: there is nothing off-diagonal.
"""
function find_offdiagonal_maximun(ih::Matrix;tol = 1E-13,zeromode=false) #use
    N = Int(size(ih)[1]/2)
    index = [1,1]
    temp_max = 0.0
    abnormal_index = [1,1]
    flag = 1

    @inbounds for i = 1:2N
        @inbounds for j = i+1:2N # iterate up-triangle part.
            t = abs(ih[i,j])
            if t>temp_max
                if zeromode && abs(ih[i,j]/(ih[j,j]+ih[i,i]))>0.4999
                    abnormal_index = [i,j]
                else
                    temp_max = t
                    index = [i,j]
                end
            end
        end
    end

    if temp_max < tol && abnormal_index==[1,1] # true: all non-zeromode elements are zero
        return 0.0,[1,1],3
    end

    if temp_max < tol && abnormal_index!=[1,1]
        index = abnormal_index
        return ih[index...],index,2
    end

    if temp_max > tol
        return ih[index...],index,1
    end
end

function given_transform!(H::Matrix,G::Matrix,i::Int,j::Int,temp_space::Matrix) # keep i<j
    N = Int(size(H)[1]/2)
    if i<=N && j>N
        given_abnormal_composite_rotation!(H,G,i,j,temp_space)
    else
        given_normal_composite_rotation!(H,G,i,j,temp_space)
    end
end

"""
This function return a prepared G to make iteration steps converge much faster
"""
function initialize_givens_eigen_solver(ih::Matrix;perturbation::Int = 1000,hamiltonian_type = "Symmetry")
    N = Int(size(ih)[1]/2)
    η = diagm(vcat([1.0 for i=1:N],[-1.0 for i=1:N]))
    # S,V = eigen(η*ih)
    r = [i/N/perturbation for i=1:N]
    r = vcat(r,r)
    # th = V*diagm(S+sort(vcat(-r,r)))*V^(-1)  #add perturbation
    th = ih + diagm(r)

    S,G = eigen(η*th)
    for i = 1:2N  # normalization
        c = transpose(G[:,i])*η*G[:,i]
        G[:,i] .=  G[:,i]./sqrt(abs(c))
    end

    for i = 1:N  # d~d^{\dagger} construction
        G[:,i] .= G[:,i+N]  #Why this construnction raise too large error?
    end
    if hamiltonian_type == "Symmetry"
        G[1:N,N+1:2N].= G[N+1:2N,1:N]
        G[N+1:2N,N+1:2N].= G[1:N,1:N]
    else hamiltonian_type == "AntiSymmetry"
        G[1:N,N+1:2N].= -G[N+1:2N,1:N]
        G[N+1:2N,N+1:2N].= G[1:N,1:N]
    end
    # the previous steps would raise small error, we need to use given eigen solver to eliminate them.
    return transpose(G)*ih*G,G
end

export deal_zeromodes

"""
    A normal output result for given_eigen_solver is ih and V (may include zero mode).
    We have two mathod to deal with it.
    type=normal:
        This function provide API to elimante zeromode by "math-well-defined rotation(physically-ill)"

            a = [1 -1; 1 1] and
            [0 0;0 4 ] = a^T*[1 -1;-1 1]*a
            [4 0;0 0 ] = a^T*[1  1; 1 1]*a

    type=abnormal:
        This function provide API to elimante zeromode by "math-ill-defined rotation(physical-well)"

            a = [1 1; 1 1] and
            [0 0;0 0] = a^T*[1 -1;-1 1]*a
            a = [1 -1; -1 1]
            [0 0;0 0] = a^T*[1  1; 1 1]*a
"""
deal_zeromodes(ih::Matrix,G::Matrix;tol=1E-13,type=:normal) = deal_zeromodes!(copy(ih),copy(G),tol=tol,type=type)

export deal_zeromodes!
function deal_zeromodes!(ih::Matrix,G::Matrix;tol=1E-13,type=:normal)
    if type == :normal
        return deal_zeromodes_normal!(ih,G,tol=tol)
    elseif type == :abnormal
        return deal_zeromodes_abnormal!(ih,G,tol=tol)
    else
        error("deal_zeromodes: Invalid type:$type")
    end
end

function deal_zeromodes_abnormal!(ih::Matrix,G::Matrix;tol=1E-13)
    for iter = 1:1000
        N = size(G)[1]
        temp_max,index,flag = find_offdiagonal_maximun(ih;tol=tol,zeromode=true)
        if flag == 3
            break
        end

        i,j = index
        sgn = sign(ih[i,i]/ih[i,j])

        tG = Matrix{Float64}(I,N,N)
        tG[i,j] = -1*sgn
        tG[j,i] = -1*sgn
        ih[i,j]=0
        ih[j,i]=0
        G .= G*tG

        if sgn<0
            ih[j,j]=0
            ih[i,i]=0
        else
            ih[i,i]=0
            ih[j,j]=0
        end
    end
    return ih,G
end

function deal_zeromodes_normal!(ih::Matrix,G::Matrix;tol=1E-13)
    for iter = 1:1000
        N = size(G)[1]
        temp_max,index,flag = find_offdiagonal_maximun(ih;tol=tol,zeromode=true)
        if flag == 3
            break
        end

        i,j = index
        sgn = sign(ih[i,i]/ih[i,j])

        tG = Matrix{Float64}(I,N,N)
        tG[i,j] = -1
        tG[j,i] = 1
        ih[i,j]=0
        ih[j,i]=0
        G .= G*tG

        if sgn<0
            ih[j,j]=4*ih[i,i]
            ih[i,i]=0
        else
            ih[i,i]=4*ih[j,j]
            ih[j,j]=0
        end
    end
    return ih,G
end

"""
    Givens Roation Eigen solver; Mathematical motivation and derivation can be found in notes dir:
"""
function given_eigen_solver(ih::Matrix ;max_iter=10000000,tol=1E-13,initialize=true,hamiltonian_type="AntiSymmetry",zeromode=false)
    N = Int(size(ih)[1]/2)
    if initialize
        ih, G = initialize_givens_eigen_solver(ih,hamiltonian_type=hamiltonian_type)
    else
        G = diagm([1 for i =1:2N])
    end

    temp_space = zeros(Float64,2N,2N)
    @inbounds for iter = 1:max_iter
        temp_max,index,flag = find_offdiagonal_maximun(ih;tol=tol,zeromode=true)
        if abs(temp_max) < tol || flag!=1
            # if zeromode && flag == 2
            #     deal_zeromodes!(ih,G) # make ih diagonalized, then G will be not canonical.
            # end
            return ih,G
        end
        i,j = index
        given_transform!(ih,G,i,j,temp_space)
    end
    print("Iteration doesnot converge!\n")
    return ih,G
end

function given_normal_composite_rotation!(H::Matrix,G::Matrix,i::Int,j::Int,temp_space::Matrix)
    N = Int(size(H)[1]/2)
    given_normal_rotation!(H,G,i,j,temp_space)
    if i > N
        given_normal_rotation!(H,G,i-N,j-N,temp_space)
    else
        given_normal_rotation!(H,G,i+N,j+N,temp_space)
    end
end

function given_abnormal_rotation!(H::Matrix,G::Matrix,i::Int,j::Int,temp_space::Matrix)
    tol_sign(x) = abs(x)>1E-11 ? sign(x) : 0
    N = Int(size(H)[1]/2)
    t = -2H[i,j]/(H[j,j]+H[i,i])
    if  abs(t) > 0.9999999 || isnan(t)
            error("Givens Rotation fail!")
            return nothing
    end
    θ = 0.5*atanh(t)
    c,s = cosh(θ),sinh(θ)
    #@show c,s
    temp_G = temp_space
    # G .= G*temp_G
    @inbounds  for id = 1:2N
        temp_space[id,i] =  c*G[id,i]
        temp_space[id,i] += s*G[id,j]
        temp_space[id,j] =  c*G[id,j]
        temp_space[id,j] += s*G[id,i]
        G[id,i] = temp_space[id,i]
        G[id,j] = temp_space[id,j]
    end
    # H .= temp_G*H*temp_G
    # #effect of temp_G
    temp_G .= H # save memry
    @inbounds  for id = 1:2N
        temp_G[i,id] = c*H[i,id]
        temp_G[i,id] += s*H[j,id]
        temp_G[j,id] = c*H[j,id]
        temp_G[j,id] += s*H[i,id]
    end
    temp_G[:,i] .= temp_G[i,:]
    temp_G[:,j] .= temp_G[j,:]
    temp_G[i,i] =  c^2*H[i,i]+2*c*s*H[i,j]+s^2*H[j,j]
    temp_G[j,j] =  c^2*H[j,j]+2*c*s*H[i,j]+s^2*H[i,i]
    temp_G[i,j] = c*s*(H[i,i]+H[j,j])+(c^2+s^2)H[i,j]
    temp_G[j,i] = temp_G[i,j]
    H .= temp_G
end

function given_abnormal_composite_rotation!(H::Matrix,G::Matrix,i::Int,j::Int,temp_space::Matrix)
    N = Int(size(H)[1]/2)
    a,b = i,j-N
    given_abnormal_rotation!(H,G,a,b+N,temp_space)
   #given_abnormal_rotation!(H,G,a+N,b,temp_space)
end

function given_normal_rotation!(H::Matrix,G::Matrix,i::Int,j::Int,temp_space::Matrix;ϵ=1E-10)
    N = Int(size(H)[1]/2)
    if abs(H[i,i]-H[j,j])<ϵ
        θ = pi/4
    else
        t = 2H[i,j]/(H[j,j]-H[i,i])
        θ = 1/2*atan(t)
    end
    c = cos(θ)
    s = sin(θ)

    temp_G = temp_space

    # G .= G*temp_G
    @inbounds  for id = 1:2N
        temp_space[id,i] =  c*G[id,i]
        temp_space[id,i] -= s*G[id,j]
        temp_space[id,j] =  c*G[id,j]
        temp_space[id,j] += s*G[id,i]
        G[id,i] = temp_space[id,i]
        G[id,j] = temp_space[id,j]
    end
    # H .= transpose(temp_G)*H*temp_G
    # @show H[i,j]
    # effect of temp_G
    temp_G .= H # save memry
    @inbounds  for id = 1:2N
        temp_G[i,id] = c*H[i,id]
        temp_G[i,id] -= s*H[j,id]
        temp_G[j,id] = c*H[j,id]
        temp_G[j,id] += s*H[i,id]
    end
    temp_G[:,i] .= temp_G[i,:]
    temp_G[:,j] .= temp_G[j,:]
    temp_G[i,i] =  c^2*H[i,i]-2*c*s*H[i,j]+s^2*H[j,j]
    temp_G[j,j] =  c^2*H[j,j]+2*c*s*H[i,j]+s^2*H[i,i]
    temp_G[i,j] = 0.0
    temp_G[j,i] = 0.0
    # display(transpose(tt)*H*tt-temp_G)
    H .= temp_G
    # H .= transpose(tt)*H*tt
end

# ReCanonicalize is not polished!!!
# include("./ReCanonicalize.jl")
# export recanonicalize

end


text:
"""
    maxwell_boltzmann(f)

Maxwell Boltzmann entropy
"""
maxwell_boltzmann(f) = f * log(f) - f

"""
    maxwell_boltzmann_prime(x)

Prim of Maxwell Boltzmann entropy
"""
maxwell_boltzmann_prime(x) = log(x)

"""
    maxwell_boltzmann_dual(f)

Dual of Maxwell Boltzmann entropy
"""
maxwell_boltzmann_dual(f) = exp(f)

"""
    maxwell_boltzmann_dual_prime(f)

Dual prim of Maxwell Boltzmann entropy
"""
maxwell_boltzmann_dual_prime(f) = exp(f)

"""
    kinetic_entropy(α, m, weights)

Reconstruct mathematical entropy from Legendre dual
"""
function kinetic_entropy(α::AbstractArray, m::AbstractArray, weights::AbstractVector)
    B = KitBase.maxwell_boltzmann_dual_prime.(α' * m)[:]
    return sum(maxwell_boltzmann.(B) .* weights)
end

ext: jl


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

mutable struct TemplateHashResult <: SwaggerModel
    minifiedTemplate::Any # spec type: Union{ Nothing, String } # spec name: minifiedTemplate
    templateHash::Any # spec type: Union{ Nothing, String } # spec name: templateHash

    function TemplateHashResult(;minifiedTemplate=nothing, templateHash=nothing)
        o = new()
        validate_property(TemplateHashResult, Symbol("minifiedTemplate"), minifiedTemplate)
        setfield!(o, Symbol("minifiedTemplate"), minifiedTemplate)
        validate_property(TemplateHashResult, Symbol("templateHash"), templateHash)
        setfield!(o, Symbol("templateHash"), templateHash)
        o
    end
end # type TemplateHashResult

const _property_map_TemplateHashResult = Dict{Symbol,Symbol}(Symbol("minifiedTemplate")=>Symbol("minifiedTemplate"), Symbol("templateHash")=>Symbol("templateHash"))
const _property_types_TemplateHashResult = Dict{Symbol,String}(Symbol("minifiedTemplate")=>"String", Symbol("templateHash")=>"String")
Base.propertynames(::Type{ TemplateHashResult }) = collect(keys(_property_map_TemplateHashResult))
Swagger.property_type(::Type{ TemplateHashResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TemplateHashResult[name]))}
Swagger.field_name(::Type{ TemplateHashResult }, property_name::Symbol) =  _property_map_TemplateHashResult[property_name]

function check_required(o::TemplateHashResult)
    true
end

function validate_property(::Type{ TemplateHashResult }, name::Symbol, val)
end


text:
# # Regression on a synthetic dataset
#
# In this page we will see how to perform Least Squares Support Vector Regression using
# `LeastSquaresSVM`. To accomplish this task, we will use synthetic data as created by the
# `make_regression` function from `MLJ`.
#
# First, we need to import all the necessary packages.
using LeastSquaresSVM
using MLJ, MLJBase
using DataFrames, CSV
using CategoricalArrays, Random
using Plots
gr();
rng = MersenneTwister(88);

# We then create the regression problem. To really push the implementation we will create
# a problem with 5 features and 500 instances/observations. We will also add a little bit
# of Gaussian noise to the problem.
X, y = MLJ.make_regression(500, 5; noise=1.0, rng=rng);

# We need to construct a `DataFrame` with the arrays created to better handle the data,
# as well as a better integration with `MLJ`.
df = DataFrame(X);
df.y = y;

# A very important part of the `MLJ` framework is its use of `scitypes`, a special kind of
# type that work together with the objects from the framework. Because the regression
# problem has the `Julia` types we need to convert this types to correct `scitypes` such
# such that the `machine`s from `MLJ` work fine.
dfnew = coerce(df, autotype(df));

# We can then observe the first three columns, together with their new types.
first(dfnew, 3) |> pretty

# We should also check out the basic statistics of the dataset.
describe(dfnew, :mean, :std, :eltype)

# Recall that we also need to standardize the dataset, we can see here that the mean is
# close to zero, but not quite, and we also need an unitary standard deviation.

# Split the dataset into training and testing sets.
y, X = unpack(dfnew, ==(:y), colname -> true);
train, test = partition(eachindex(y), 0.75, shuffle=true, rng=rng);
stand1 = Standardizer();
X = MLJBase.transform(MLJBase.fit!(MLJBase.machine(stand1, X)), X);

# We should make sure that the features have mean close to zero and an unitary standard
# deviation.
describe(X |> DataFrame, :mean, :std, :eltype)

# Define a good set of hyperparameters for this problem and train the regressor. We will
# use the amazing capability of `MLJ` to tune a model and return the best model found.
#
# For this we have taken some judicious guessing on the best values that the hyperparameters
# should take. We employ 5-fold cross-validation and a 400 by 400 grid of points to do
# an exhaustive search.
#
# We will train the regressor using the root mean square error which is defined as follows
#
# ```math
# RMSE = \sqrt{\frac{\sum_{i=1}^{N} \left(\hat{y}_i - y_i \right)^2}{N}}
# ```
#
# where we define $\hat{y}_i$ as the *predicted value*, and $y_i$ as the real value.
model = LeastSquaresSVM.LSSVRegressor();
r1 = MLJBase.range(model, :σ, lower=7e-4, upper=1e-3);
r2 = MLJBase.range(model, :γ, lower=120.0, upper=130.0);
self_tuning_model = TunedModel(
    model=model,
    tuning=Grid(goal=400, rng=rng),
    resampling=CV(nfolds=5),
    range=[r1, r2],
    measure=MLJBase.rms,
    acceleration=CPUThreads(), # We use this to enable multithreading
);

# And now we proceed to train all the models and find the best one!
mach = MLJ.machine(self_tuning_model, X, y);
MLJBase.fit!(mach, rows=train, verbosity=0);
fitted_params(mach).best_model

# Having found the best hyperparameters for the regressor model we proceed to check how the
# model generalizes and we use the test set to check the performance.
ŷ = MLJBase.predict(mach, rows=test);
result = round(MLJBase.rms(ŷ, y[test]), sigdigits=4);
@show result # Check the result

# We can see that we did quite well. A value of 1, or close enough, is good. We expect it
# to reach a lower value, closer to zero, but maybe we needed more refinement in the grid
# search.
#
# We can also see a plot of the predicted and true values. The closer these dots are to the
# diagonal means that the model performed well.
scatter(ŷ, y[test], markersize=9)
r = range(minimum(y[test]), maximum(y[test]); length=length(test))
plot!(r, r, linewidth=9)
plot!(size=(3000, 2100))

# We can actually see that we are not that far off, maybe a little more search could
# definitely improve the performance of our model.

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
# Autogenerated wrapper script for Usrsctp_jll for x86_64-linux-gnu
export libusrsctp

JLLWrappers.@generate_wrapper_header("Usrsctp")
JLLWrappers.@declare_library_product(libusrsctp, "libusrsctp.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libusrsctp,
        "lib/libusrsctp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@inline function find_project_file(package_path)
    for project_file in Base.project_names
        project_file_path = joinpath(package_path, project_file)
        pkg = Pkg.Types.read_project(project_file_path)
        if pkg.name !== nothing
            return project_file_path
        end
    end
    throw(ErrorException("Could not find project file in package at \"$(package_path)\""))
end


text:
module JSONSchema2Types

using JSON3

export generate_type_module

include("parser.jl")

end


text:
include("reactor.jl")
include("scan.jl")
include("study.jl")
include("sampling.jl")
include("equation.jl")
include("equation_set.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "GStreamer"
version = v"1.18.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://gstreamer.freedesktop.org/src/gstreamer/gstreamer-$(version).tar.xz", "0c2e09e18f2df69a99b5cb3bd53c597b3cc2e35cf6c98043bb86a66f3d312100")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
apk add bash-completion libcap gettext
cd gstreamer-*
mkdir build
cd build
meson .. --cross-file=${MESON_TARGET_TOOLCHAIN}
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(exclude=[Platform("i686", "linux", libc = "musl"), Platform("powerpc64le", "linux")])

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgstbase-1.0", "libgstbase-1"], :libgstbase),
    LibraryProduct(["libgstcheck-1.0", "libgstcheck-1"], :libgstcheck),
    LibraryProduct(["libgstcontroller-1.0", "libgstcontroller-1"], :libgstcontroller),
    LibraryProduct(["libgstnet-1.0", "libgstnet-1"], :libgstnet),
    LibraryProduct(["libgstreamer-1.0", "libgstreamer-1"], :libgstreamer),
    ExecutableProduct("gst-inspect-1.0", :gst_inspect),
    ExecutableProduct("gst-launch-1.0", :gst_launch),
    ExecutableProduct("gst-stats-1.0", :gst_stats),
    ExecutableProduct("gst-tester-1.0", :gst_tester),
    ExecutableProduct("gst-typefind-1.0", :gst_typefind),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Glib_jll", uuid="7746bdde-850d-59dc-9ae8-88ece973131d"), v"2.59.0"; compat="2.59.0")
    Dependency(PackageSpec(name="LibUnwind_jll", uuid="745a5e78-f969-53e9-954f-d19f2f74f4e3"))
    Dependency(PackageSpec(name="Elfutils_jll", uuid="ab5a07f8-06af-567f-a878-e8bb879eba5a"))
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"), v"6.1.2")
    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
    Dependency(PackageSpec(name="libcap_jll", uuid="eef66a8b-8d7a-5724-a8d2-7c31ae1e29ed"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# # Maximum Likelihood Estimation

using Distributions
using HMMBase
using Plots

y1 = rand(Normal(0,2), 1000)
y2 = rand(Normal(10,1), 500)
y = vcat(y1, y2, y1, y2)

plot(y, linetype=:steppre, size=(600,200))

# For now `HMMBase` does not handle the initialization of the parameters.
# Hence we must instantiate an initial HMM by hand.

hmm = HMM([0.5 0.5; 0.5 0.5], [Normal(-1,2), Normal(15,2)])
hmm = fit_mle(hmm, y, display = :iter)
#-

z_viterbi = viterbi(hmm, y)
plot(z_viterbi, linetype=:steppre, label="Viterbi decoded hidden state", size=(600,200))


text:
"""
    sample_tips(tree::Tree, f::Float64 ; root_age = true)

Sample tips from a tree. Each tip is kept with probability `f`. The fuction returns a tupple `(phylo, rates)`, where phylo is a `Tree` object containing the sampled tree, and rates a vector of `Floats` with its tip rates.

# Arguments
- `tree::Tree`: the phylogeny.
- `f::Float64`: the sampling probability.

# Keyword arguments
- `root_age::Bool`: If true, the sampling is constrained so that the tree keeps the same root age : at least one tip is sampled in each of the two tree subtrees.
"""
function sample_tips(tree::Tree, f::Float64 ; root_age = true)
    function combine(left, right, bl, att)
        extant = (left.extant || right.extant)
        return Tree([left,right], bl, att, extant, "")
    end
    function aux(subtree)
        u = 0.5
        if subtree.extant
            if length(subtree.offsprings) == 0
                if f < 1.
                    u = rand()
                end
                if u < f
                    return subtree
                else
                    return Tree()
                end
            else
                combine(aux(subtree.offsprings[1]), aux(subtree.offsprings[2]), subtree.branch_length, subtree.attributes)
            end
        else
            return subtree
        end
    end

    if root_age
        left_tree = Tree()
        while (!left_tree.extant)
            left_tree = aux(tree.offsprings[1])
        end
        right_tree = Tree()
        while (!right_tree.extant)
            right_tree = aux(tree.offsprings[2])
        end
        new_tree = Tree([left_tree,right_tree], tree.branch_length, tree.attributes, tree.label)
        tip_rates = extract_tip_rates_light(new_tree, return_extinct = false)
        return (prune_extinct_lineages(new_tree), tip_rates)
    else
        new_tree = aux(tree)
        tip_rates = extract_tip_rates_light(new_tree, return_extinct = false)
        return (prune_extinct_lineages(new_tree), tip_rates)
    end
end

function extract_tip_rates_light(tree::Tree ; id = 1, return_extinct = true)

    function aux(subtree, x)
        if subtree.n_nodes == 0
        elseif length(subtree.offsprings) == 0
            if subtree.extant || return_extinct
                pushfirst!(x,subtree.attributes[id])
            end
        else
            aux(subtree.offsprings[2], x)
            aux(subtree.offsprings[1], x)
        end
    end

    x = Array{Float64,1}(undef,0)
    aux(tree,x)
    return x
end

function prune_extinct_lineages(tree::Tree)
    function combine(left, right, bl, att)
        if ! left.extant
            if ! right.extant
                return Tree()
            else
                return Tree(right.offsprings, right.branch_length +bl,
                    att, true, right.label)
            end
        elseif ! right.extant
            return Tree(left.offsprings, left.branch_length + bl,
                att, true, left.label)
        else
            return Tree([left, right], bl, att, true, "")
        end
    end

    function aux(subtree)
        if length(subtree.offsprings) == 0
            return subtree
        else
            combine(aux(subtree.offsprings[1]), aux(subtree.offsprings[2]), subtree.branch_length, subtree.attributes)
        end
    end

    aux(tree)
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
push!(LOAD_PATH, "Modules")
using Colors
using JuliRay
Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))

function Smooth(x::Real)
    f(x) = ifelse(x>0,exp(-1/x),0)
    return f(2x/√3)/(f(2x/√3)+f(2/√3-2x/√3))
end
function Smooth(a::Real, b::Real, x::Real)
    return Smooth((x-a)/(b-a))
end
function Smooth(dict::Dict{A,B} where A <: Real where B <: Any, t::Real)
    T=sort(collect(keys(dict)))
    L=length(T)
    V=[dict[T[i]] for i in 1:L]

    if t ≤ T[1]
        return V[1]
    end
    for i in 2:L
        if t ≤ T[i]
            return V[i-1]+(V[i]-V[i-1])*Smooth(T[i-1],T[i],t)
        end
    end
    return V[end]
end

n=5
colors=[
    RGB(0.8,0.05,0.15), # Red
    RGB(0.8,0.3,0), # Orange
    RGB(0.8,0.7,0), # Yellow
    RGB(0,0.8,0.3), # Mint Green
    RGB(0,0.6,0.8)  # Cyan Blue
    ]

# S(i)=rgbColor(Box([0,0,n-i].+0.03,[i,i,n-i+1].-0.03), colors[i])
S(i)=rgbColor(RoundedBox([0,0,n-i].+0.001,[i,i,n-i+1].-0.001, 0.1), colors[i])

Ss1(t)=csgUnion(S(i) for i in 1:n)
Ss2(t)=ParallelTranslation(Rotate(Ss1(t),[1,1,-1], 2*π/3, fixedpoint=[1,1,n]), [0,Smooth(Dict(0=>2n, 1=>0),t),Smooth(Dict(0=>1.5n, 1=>-1),t)])
Ss3(t)=ParallelTranslation(Rotate(Ss1(t),[1,1,-1], -2*π/3, fixedpoint=[1,1,n]), [-1,0,Smooth(Dict(1=>2n,2=>0),t)])
SSS(t)=ParallelTranslation(
    Rotate(
        csgUnion(Ss1(t), Ss2(t), Ss3(t)),[0,1,0], Smooth(Dict(2=>0°, 4=>-90°),t), fixedpoint=[n,n,n]/2
    ),  Smooth(Dict(2=>[0,0,0], 4=>[4,-4,0], 5=>[-3,-3,0]),t)
)

SSS′(t)=Rotate(
    ParallelTranslation(
        Rotate(
            ParallelTranslation(
                Scaling(SSS(4),-1), [0,0,n]
            ), [0,0,1], -90°
        ), [-1,1,0]*Smooth(Dict(2=>3n, 4=>0),t)+[-1,0,0]
    ), [0,0,1], Smooth(Dict(4=>0, 5=>90°),t)
)

N = 20
for i ∈ 1:5N
    t=i/N
    lng=Smooth(Dict(0=>60°, 1=>30°, 2=>-20°, 4=>-150°, 5=>-120°),t)
    lat=20°
    lookat=[n,n,n]/2+Smooth(Dict(4=>[0,0,0], 5=>[-5.5,0.5-2,0]),t)
    zoom=Smooth(Dict(2=>0.1, 4=>0.055, 5=>0.07),t)
    render(csgUnion(SSS(t),SSS′(t)),camera=LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=zoom,lookat=lookat),name="SumOfSquare",index=i)
end

for i ∈ 5N:7N
    t=i/N
    s1=ParallelTranslation(Rotate(Ss1(0), [0,0,1], 180°),[-3,3,0])
    s2=ParallelTranslation(Rotate(s1,[1,1,1],120°, fixedpoint=[-4,2,5]),[0,0,-1])
    s3=ParallelTranslation(Rotate(s1,[1,1,1],-120°, fixedpoint=[-4,2,5]),[-1,0,-1])
    s4=Scaling(s1,-1,fixedpoint=[-3.5, 0, 2.5])
    s5=Scaling(s2,-1,fixedpoint=[-3.5, 0, 2.5])
    s6=Scaling(s3,-1,fixedpoint=[-3.5, 0, 2.5])

    t1=s1
    t2=Transparent(s2,FT(0,Smooth(Dict(6.0=>0, 6.5=>1),t)))
    t3=Transparent(s3,FT(0,Smooth(Dict(5.0=>0, 5.5=>1),t)))
    t4=Transparent(s4,FT(0,Smooth(Dict(5.5=>0, 6.0=>1),t)))
    t5=Transparent(s5,FT(0,Smooth(Dict(6.0=>0, 6.5=>1),t)))
    t6=Transparent(s6,FT(0,Smooth(Dict(6.5=>0, 7.0=>1),t)))
    obj=csgUnion(t1,t2,t3,t4,t5,t6)
    lng=Smooth(Dict(0=>60°, 1=>30°, 2=>-20°, 4=>-150°, 5=>-120°, 7=>-120°),t)
    lat=20°
    lookat=[n,n,n]/2+Smooth(Dict(4=>[0,0,0], 5=>[-5.5,0.5-2,0], 7=>[-8,-2,0]),t)
    zoom=Smooth(Dict(2=>0.1, 4=>0.055, 5=>0.07, 7=>0.1),t)
    render(obj,camera=LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=zoom,lookat=lookat),name="SumOfSquare",index=i)
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
@doc """
""" ->
function shuffle(data)
    m,n = size(data)
    inds = randperm(m)
    data[inds,:]
end

function split_class_data(data)
    m,n = size(data)
    point1 = round(Int64, 0.6*m)
    point2 = round(Int64, 0.8*m)
    train = data[1:point1,:]
    val   = data[point1+1:point2,:]
    test  = data[point2+1:end,:]

    train,val,test
end

@doc """ split each class dataset
""" ->
function split_data(data)
    num_label = maximum(data[:,end])

    function reducer(data1, data2)
        (vcat(data1[1],data2[1]), vcat(data1[2], data2[2]), vcat(data1[3], data2[3]))
    end

    data = @parallel (reducer) for i = 1:num_label
        split_class_data(data[data[:,end] .== i,:])
    end

    data
end

@doc """ eval accuracy of model
""" ->
function acc(preds,gts)
    @assert length(preds) == length(gts)
    sum(preds .== gts) / length(gts)
end

@doc """ xgboost for train,val,test dataset without gridtune
         eval with acc.
""" ->
function model_xgboost(train, val, test, num_iter::Int64)
    train_X,train_Y = train[:,1:end-1],vec(full(train[:,end]))
    val_X,  val_Y = val[:,1:end-1],vec(full(val[:,end]))
    test_X, test_Y = test[:,1:end-1],vec(full(test[:,end]))

    num_class = maximum(train_Y)+1
    m_tr,n_tr   = size(train_X)
    m_val,n_val = size(val_X)
    @info("in the model, there are $num_class classes")
    @info("train_X size is :$m_tr,$n_tr")
    @info("val_X size is : $m_val, $n_tr")
    #stat_class_sample(train,class_names())
    #sleep(3)

    if debug
        num_round = 1
    else
        num_round = num_iter
    end

    dtrain = DMatrix(train_X, label = train_Y)
    dval   = DMatrix(val_X,   label = val_Y)
    dtest  = DMatrix(test_X,  label = test_Y)

    watch_list = [(dtrain,"train"), (dval,"val")]
    param      = Dict{ASCIIString,Any}("max_depth"=>2,"eta"=>0.01,"nthread"=>68,
                                       "objective"=>"multi:softmax","silent"=>1,
                                       "alpha"=>0.85,
                                       "sub_sample"=>0.6,"num_class"=>num_class)

    bst = xgboost(dtrain,num_round,watchlist=watch_list,param=param,metrics=["merror"],seed=2015)
    save(bst,"xgb_night.model")
    test_preds = XGBoost.predict(bst, test_X)
    mean_f1 = acc(test_preds, get_info(dtest,"label"))
    @info("on test dataset mean-f1-score is $mean_f1")
    test_preds = map(x->convert(Int64,x), test_preds)
    test_Y = map(x->convert(Int64,x), test_Y)
    num_class = convert(Int64, num_class)

    ROC_print(test_Y.+1,test_preds.+1,num_class = num_class, challenge_list=class_names())

end


text:
# Utils
function module_dir()
    return joinpath(@__DIR__, "..")
end

# Activate package
using Pkg
Pkg.activate(module_dir())

using LinearAlgebra
using LyceumBase
using LyceumMuJoCo
using LyceumMuJoCoViz
using MuJoCo
using Plots
using JLD2
using Random

include("methods.jl")

################################################################################
# Demo
################################################################################

timestep=0.01
jm, jd, mjsim = mj_model("box.xml", timestep=timestep)
tsim = 3.0
N = Int(floor(tsim/timestep))

jd.qpos .= [-1.5,-0.50,0.50,1,0,0,0]
jd.qvel .= [4.0,0.8,0,0,0,0]

traj = zeros(20, 0)
tcompute = 0.0
ztraj = []
for i = 1:N
    jd.ctrl .= 0.0
    tcompute += @elapsed mj_step(jm, jd);
    traj = hcat(traj, deepcopy(getstate(mjsim)))
	x = deepcopy(jd.qpos[1:3])
	v = deepcopy(jd.qvel[1:3])
	q = deepcopy(jd.qpos[4:7])
	ω = deepcopy(jd.qvel[4:6])
	push!(ztraj, [x; v; q; ω])
end
tsim / tcompute

LyceumMuJoCoViz.visualize(mjsim, trajectories=traj)

jd.qpos
jd.qvel

# Saving trajectory
path = joinpath(module_dir(), "results/cone_compare.jld2")
jldsave(path, ztraj=ztraj)


text:
module JSDeps
using JSON

struct JSRemoteLibrary
    name :: String
    url :: String
    # The following become part of the shim
    deps :: Union{Nothing, Vector{String}}
    exports :: Union{Nothing, String}
end

function shimdict(lib::JSRemoteLibrary)
    isnothing(lib.deps) && isnothing(lib.exports) && return nothing
    shim = Dict{Symbol,Any}()
    if !isnothing(lib.deps)
        shim[:deps] = lib.deps
    end
    if !isnothing(lib.exports)
        shim[:exports] = lib.exports
    end
    return shim
end

struct JSSnippet
    deps :: Vector{String}
    args :: Vector{String}
    js :: String
end

macro js_str(s)
    return s
end

struct RequireJS
    libraries :: Dict{String, JSRemoteLibrary}
    snippets :: Vector{JSSnippet}

    RequireJS() = new(Dict(), [])
end

function shimdict(r::RequireJS)
    shim = Dict{String,Any}()
    for (name, lib) in r.libraries
        @assert name == lib.name
        libshim = shimdict(lib)
        if !isnothing(libshim)
            shim[name] = libshim
        end
    end
    return shim
end

function Base.push!(r::RequireJS, lib::JSRemoteLibrary)
    if lib.name in keys(r.libraries)
        error("Library already added.")
    end
    r.libraries[lib.name] = lib
end

Base.push!(r::RequireJS, s::JSSnippet) = push!(r.snippets, s)

function verify(r::RequireJS)
    for s in r.snippets
        for dep in s.deps
            dep in keys(r.libraries) || error("$(dep) missing from libraries")
        end
    end
end

function writejs(filename::AbstractString, r::RequireJS)
    open(filename, "w") do io
        writejs(io, r)
    end
end

function writejs(io::IO, r::RequireJS)
    write(io, """
    // Generated by Documenter.jl
    requirejs.config({
      paths: {
    """)
    for (name, lib) in r.libraries
        url = endswith(lib.url, ".js") ? replace(lib.url, r"\.js$" => "") : lib.url
        write(io, """
            '$(lib.name)': '$(url)',
        """) # FIXME: escape bad characters
    end
    write(io, """
      },
    """)

    shim = shimdict(r)
    if !isempty(shim)
        write(io, "  shim: ")
        JSON.print(io, shim, 2) # FIXME: escape JS properly
        write(io, ",\n")
    end

    write(io, """
    });
    """)

    for s in r.snippets
        args = join(s.args, ", ") # FIXME: escapes
        deps = join(("\'$(d)\'" for d in s.deps), ", ") # FIXME: escapes
        write(io, """
        require([$(deps)], function($(args)) {
        $(s.js)
        })
        """)
    end
end

"""
    read_js_snippet(filename::AbstractString) -> JSSnippet
    read_js_snippet(io::IO) -> JSSnippet

Parses a JS snippet file into a [`JSSnippet`](@ref) object.
"""
function read_js_snippet end

read_js_snippet(filename::AbstractString; kwargs...) = open(filename, "r") do io
    read_js_snippet(io; kwargs...)
end

function read_js_snippet(io::IO)
    libraries = String[]
    arguments = String[]
    lineno = 1
    while true
        pos = position(io)
        line = readline(io)
        m = match(r"^//\s*([a-z]+):(.*)$", line)
        if m === nothing
            seek(io, pos) # undo the last readline() call
            break
        end
        if m[1] == "libraries"
            libraries = strip.(split(m[2], ","))
            if any(s -> match(r"^[a-z-_]+$", s) === nothing, libraries)
                error("Unable to parse a library declaration '$(line)' on line $(lineno)")
            end
        elseif m[1] == "arguments"
            arguments = strip.(split(m[2], ","))
        end
        lineno += 1
    end
    snippet = String(read(io))
    JSAPI.JSSnippet(libraries, arguments, snippet)
end

# struct JSSnippetParseError
#     key :: String
#     line :: String
#     lineno :: Integer
# end
#
# function Base.showerror(io::IO, err::JSSnippetParseError)
#     println(io, "JSSnippetParseError: bad '$(err.key)' declaration on line $(err.lineno)")
#     println(io, "  attempted to parse: '$(err.line)'")
# end

end


text:
export CuGraph, capture, instantiate, CuGraphExec, launch, update,
       capture_status, is_capturing,
       @captured

## graph

@enum_without_prefix CUstreamCaptureMode CU_

mutable struct CuGraph
    handle::CUgraph
    ctx::CuContext

    function CuGraph(flags=STREAM_CAPTURE_MODE_GLOBAL)
        handle_ref = Ref{CUgraph}()
        cuGraphCreate(handle_ref, flags)

        ctx = CuCurrentContext()
        obj = new(handle_ref[], ctx)
        finalizer(unsafe_destroy!, obj)
        return obj
    end

    """
        capture([flags], [throw_error::Bool=true]) do
            ...
        end

    Capture a graph of CUDA operations. The returned graph can then be instantiated and
    executed repeatedly for improved performance.

    Note that many operations, like initial kernel compilation or memory allocations,
    cannot be captured. To work around this, you can set the `throw_error` keyword to false,
    which will cause this function to return `nothing` if such a failure happens. You can
    then try to evaluate the function in a regular way, and re-record afterwards.

    See also: [`instantiate`](@ref).
    """
    global function capture(f::Function; flags=STREAM_CAPTURE_MODE_GLOBAL, throw_error::Bool=true)
        cuStreamBeginCapture_v2(stream(), flags)

        ctx = CuCurrentContext()
        obj = nothing
        try
            f()
        finally
            handle_ref = Ref{CUgraph}()
            err = unsafe_cuStreamEndCapture(stream(), handle_ref)
            if err == ERROR_STREAM_CAPTURE_INVALIDATED && !throw_error
                return nothing
            elseif err != CUDA_SUCCESS
                throw_api_error(err)
            end

            obj = new(handle_ref[], ctx)
            finalizer(unsafe_destroy!, obj)
        end
        return obj::CuGraph
    end
end

function unsafe_destroy!(graph::CuGraph)
    @finalize_in_ctx graph.ctx cuGraphDestroy(graph)
end

Base.unsafe_convert(::Type{CUgraph}, graph::CuGraph) = graph.handle

## instantiated graph

mutable struct CuGraphExec
    handle::CUgraphExec
    graph::CuGraph
    ctx::CuContext

    """
        instantiate(graph::CuGraph)

    Creates an executable graph from a graph. This graph can then be launched, or updated
    with an other graph.

    See also: [`launch`](@ref), [`update`](@ref).
    """
    global function instantiate(graph::CuGraph)
        handle_ref = Ref{CUgraphExec}()
        error_node = Ref{CUgraphNode}()
        buflen = 256
        buf = Vector{UInt8}(undef, buflen)

        GC.@preserve buf begin
            cuGraphInstantiate_v2(handle_ref, graph, error_node, pointer(buf), buflen)
            diag = String(buf)
            # TODO: how to use these?
        end

        ctx = CuCurrentContext()
        obj = new(handle_ref[], graph, ctx)
        finalizer(unsafe_destroy!, obj)
        return obj
    end
end

function unsafe_destroy!(exec::CuGraphExec)
    @finalize_in_ctx exec.ctx cuGraphDestroy(exec)
end

Base.unsafe_convert(::Type{CUgraphExec}, exec::CuGraphExec) = exec.handle

"""
    launch(exec::CuGraphExec, [stream::CuStream])

Launches an executable graph, by default in the currently-active stream.
"""
launch(exec::CuGraphExec, stream::CuStream=stream()) = cuGraphLaunch(exec, stream)

@enum_without_prefix CUgraphExecUpdateResult CU_

"""
    update(exec::CuGraphExec, graph::CuGraph; [throw_error::Bool=true])

Check whether an executable graph can be updated with a graph and perform the update if
possible. Returns a boolean indicating whether the update was successful. Unless
`throw_error` is set to false, also throws an error if the update failed.
"""
function update(exec::CuGraphExec, graph::CuGraph; throw_error::Bool=true)
    error_node = Ref{CUgraphNode}()
    update_result = Ref{CUgraphExecUpdateResult}()
    cuGraphExecUpdate(exec, graph, error_node, update_result)
    if update_result[] != GRAPH_EXEC_UPDATE_SUCCESS && !throw_error
        return false
    elseif update_result[] != GRAPH_EXEC_UPDATE_SUCCESS
        error("Could not update the executable graph: $(update_result[])")
    end
    return true
end

## global properties

@enum_without_prefix CUstreamCaptureStatus CU_

function capture_status(stream::CuStream=stream())
    status_ref = Ref{CUstreamCaptureStatus}()
    id_ref = Ref{UInt64}()
    cuStreamGetCaptureInfo(stream, status_ref, id_ref)
    return (status=status_ref[],
            id=(status_ref[] == STREAM_CAPTURE_STATUS_ACTIVE ? id_ref[] : nothing))
end

is_capturing(stream::CuStream=stream()) =
    capture_status(stream).status != STREAM_CAPTURE_STATUS_NONE

## convenience macro

"""
    for ...
        @capture begin
            # code that executes several kernels or CUDA operations
        end
    end

A convenience macro for recording a graph of CUDA operations and automatically cache and
update the execution. This can improve performance when executing kernels in a loop, where
the launch overhead might dominate the execution.

!!! warning

    For this to be effective, the kernels and operations executed inside of the
    captured region should not signficantly change across iterations of the loop. It is
    allowed to, e.g., change kernel arguments or inputs to operations, as this will be
    processed by updating the cached executable graph. However, significant changes will
    result in an instantiation of the graph from scratch, which is an expensive operation.

See also: [`capture`](@ref).
"""
macro captured(ex)
    @gensym exec
    @eval __module__ begin
        const $exec = Ref{CuGraphExec}()
    end
    quote
        executed = false

        # capture
        graph = capture(throw_error=false) do
            $(esc(ex))
        end
        if graph === nothing
            # if the capture failed, this may have been due to JIT compilation.
            # execute the body out of capture, and try capturing again.
            $(esc(ex))
            graph = capture(throw_error=true) do
                # don't tolerate capture failures now so that the user will be informed
                $(esc(ex))
            end
            executed = true
        end

        # update or instantiate
        if !isassigned($(esc(exec))) || !update($(esc(exec))[], graph; throw_error=false)
            $(esc(exec))[] = instantiate(graph)
        end

        # execute
        executed || launch($(esc(exec))[])
    end
end


text:
function KS_solve_MGC_cg( pw::PWGrid, V_ionic, Focc, Nstates::Int;
                           psi0=nothing, Potentials0 = nothing,
                           α_t = 3e-5, NiterMax=1000, verbose=false )

    Ns = pw.Ns
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx
    Ω = pw.Ω

    if psi0 == nothing
        srand(2222)
        psi = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)
        psi = ortho_gram_schmidt(psi)
    else
        psi = copy(psi0)
    end

    Potentials = PotentialsT( V_ionic, zeros(Npoints), zeros(Npoints) )

    d = zeros(ComplexF64, Ngwx, Nstates)
    g_old = zeros(ComplexF64, Ngwx, Nstates)
    d_old = zeros(ComplexF64, Ngwx, Nstates)
    Kg = zeros(ComplexF64, Ngwx, Nstates)
    Kg_old = zeros(ComplexF64, Ngwx, Nstates)

    β        = 0.0
    Etot     = calc_E_MGC( pw, Potentials, psi )
    Etot_old = Etot

    ctmp = zeros(ComplexF64, Npoints, Nstates)
    idx = pw.gvecw.idx_gw2r

    S = zeros(ComplexF64,Nstates,Nstates)
    Q = zeros(ComplexF64,Nstates,Nstates)

    rho = zeros(Float64,Npoints)

    for iter = 1:NiterMax

        g = calc_grad_MGC( pw, Potentials, Focc, psi)
        Kg = Kprec(pw,g)

        println("sum(g) = ", sum(g))
        println("sum(Kg) = ", sum(Kg))

        if iter != 1
            #β = real(sum(conj(g).*Kg))/real(sum(conj(g_old).*Kg_old))
            β = real(sum(conj(g-g_old).*Kg))/real(sum(conj(g_old).*Kg_old))
            #β = real(sum(conj(g-g_old).*Kg))/real(sum(conj(g-g_old).*d))
            #β = real(sum(conj(g).*Kg))/real(sum((g-g_old).*conj(d_old)))
        end

        d = -Kprec(pw, g) + β * d_old

        Δ = 2.0*real( sum( conj(d).*g ) )
        psi_t = ortho_gram_schmidt( psi + α_t*d )
        E_trial = calc_E_MGC( pw, Potentials, psi_t )
        curvature = ( E_trial - ( Etot + α_t*Δ ) ) /α_t^2
        α = -Δ/(2*curvature)

        println("β = ", β)
        println("α = ", α)

        # Update wavefunction
        psi = psi[:,:] + α*d[:,:]

        Etot = calc_E_MGC( pw, Potentials, psi )

        diff = abs(Etot-Etot_old)
        @printf("CG step %8d = %18.10f %18.10f\n", iter, Etot, diff)
        if diff < 1e-6
            @printf("CONVERGENCE ACHIEVED\n")
            break
        end

        g_old = copy(g)
        d_old = copy(d)
        Kg_old = copy(Kg)
        Etot_old = Etot
    end
    return psi, Energies, Potentials
    #
end

function calc_E_MGC( pw::PWGrid, Potentials, ϕ::Array{ComplexF64,2}; Nel=1.0, η=1.0 )

    Ngwx    = size(ϕ)[1]
    Nstates = size(ϕ)[2]
    Npoints = prod(pw.Ns)
    Ω = pw.Ω

    S = zeros(ComplexF64,Nstates,Nstates)
    Q = zeros(ComplexF64,Nstates,Nstates)

    for j = 1:Nstates
    for i = 1:Nstates
        S[i,j] = sum( conj(ϕ[:,i]).*ϕ[:,j] )
        if i == j
            Q[i,j] = 2 - S[i,j]
        end
    end
    end

    E_kin = 0.0
    for j = 1:Nstates
    for i = 1:Nstates
        E_kin = E_kin + real(Q[i,j]*sum(conj(ϕ[:,j] .* op_K(pw, ϕ[:,i]))))
    end
    end

    ctmp = zeros(ComplexF64, Npoints, Nstates)
    idx = pw.gvecw.idx_gw2r
    for ic = 1:Nstates
        ctmp[idx,ic] = ϕ[:,ic]
    end
    ϕ_r = G_to_R( pw.Ns, ctmp )
    # orthonormalization in real space
    ortho_gram_schmidt!( Nstates, ϕ_r )
    scale!( sqrt(Npoints/Ω), ϕ_r )

    ρ = zeros(Float64,Npoints)

    for j = 1:Nstates
    for i = 1:Nstates
        for ip = 1:Npoints
            # XXX need to use conj??? Use factor 2 ??
            ρ[ip] = real( Q[i,j] * conj(ϕ_r[ip,i]) * ϕ_r[ip,j] )
        end
    end
    end

    V_Hartree = real( G_to_R( pw.Ns, Poisson_solve(pw, ρ) ) )
    E_Hartree = 0.5*dot( V_Hartree, ρ ) * Ω/Npoints

    V_XC = excVWN( ρ ) + ρ .* excpVWN( ρ )

    E_xc = dot( excVWN(ρ), ρ ) * Ω/Npoints

    E_Ionic = dot( Potentials.Ionic, ρ ) * Ω/Npoints

    intRho = sum(ρ)*Ω/Npoints
    #println("intRho = ", intRho)

    E_MGC = E_kin + E_Hartree + E_Ionic + E_xc + η*(Nel - intRho)

    #println("E_MGC Kin     = ", E_kin)
    #println("E_MGC Ionic   = ", E_Ionic)
    #println("E_MGC Hartree = ", E_Hartree)
    #println("E_MGC XC      = ", E_xc)
    #println("E_MGC Total   = ", E_MGC)

    return E_MGC

end

function calc_grad_MGC( pw::PWGrid, Potentials, Focc, psi::Array{ComplexF64,2} )

    # Focc are assumed to be 2s
    Ngwx    = size(psi)[1]
    Nstates = size(psi)[2]

    g = zeros(ComplexF64,Ngwx,Nstates)
    H = zeros(ComplexF64,Nstates,Nstates)
    S = zeros(ComplexF64,Nstates,Nstates)
    Q = zeros(ComplexF64,Nstates,Nstates)

    Hpsi = op_H( pw, Potentials, psi )

    for m = 1:Nstates
    for n = 1:Nstates
        H[m,n] = sum( conj(psi[:,m]).*Hpsi[:,n] )
        S[m,n] = sum( conj(psi[:,m]).*psi[:,n] )
        if m == n
            Q[m,n] = 2 - S[m,m]
        end
    end
    end
    #
    η = 1.0  # ????
    #
    for n = 1:Nstates
        for m = 1:Nstates
            g[:,n] = Hpsi[:,m]*Q[m,n] - psi[:,m]*Q[m,n]*η - psi[:,m]*(H[m,n] - η*S[m,n])
        end
        g[:,n] = 2*g[:,n]
    end
    #
    return g
end


text:
import Statistics: mean
using DataStructures

"""
Container for generated Mesh
"""
mutable struct Mesh
    nodes :: Vector{Vector{Real}}
    elements :: Vector{Vector{Integer}}
    elementsets :: DefaultDict{String, Vector{Integer}}
end

function emptyMesh()
    Mesh(Vector{Float64}[], Vector{Int64}[], DefaultDict{String, Vector{Integer}}(() -> Vector{Integer}()))
end

function translateNodes!(mesh::Mesh, node_ids::Vector{Int64}, translation::Vector{Float64}; copynodes::Bool=true)
    if !copynodes
        throw(DomainError())
    end

    coords = Vector{Vector{Float64}}()

    for node_id in node_ids
        push!(coords, mesh.nodes[node_id]+translation)
    end

    node_ids_new = appendNodes!(mesh, coords)

    return node_ids_new
end

function appendNodes!(mesh::Mesh, nodecoords::Vector{Vector{Float64}})
    start_id = size(mesh.nodes, 1)+1
    length = size(nodecoords, 1)

    append!(mesh.nodes, nodecoords)

    return collect(start_id:1:start_id+length-1)
end

function appendNodes!(mesh::Mesh, nodecoords::Array{Float64,2})
    start_id = size(mesh.nodes, 1)+1
    length = size(nodecoords, 1)

    for i in range(1, stop=length, step=1)
        push!(mesh.nodes, nodecoords[i,:])
    end

    return collect(start_id:1:start_id+length-1)
end

function appendElement!(mesh::Mesh, node_ids, elem_type)
    el_id = size(push!(mesh.elements, node_ids), 1)
    push!(mesh.elementsets[elem_type], el_id)
    el_id
end

function getFac(boundary::AbstractBoundary; flipdirection::Bool=false)

    N = boundary.node_num-1

    if abs(boundary.bias) != 1.0
        α = abs(boundary.bias)

        us = [[(α.^i-1)/(α^N-1) for i in range(0, stop=N, step=1)]...]

    else
        us = collect(0:1:N)/(N)
    end

    if (boundary.bias > 0) != (!flipdirection)
        return us
    else
        return (1 .- us)[end:-1:1]
    end
end

function meshBoundary!(mesh::Mesh, boundary::StraightBoundary, vertices, skipFirst::Bool, skipLast::Bool)

    # get coordinates of start and end point
    coords1 = vertices[boundary.vertice_start]
    coords2 = vertices[boundary.vertice_end]

    # Stretching functions
    Δcoords = coords2 - coords1

    # distribution
    us = getFac(boundary)
    nodes = coords1' .+ Δcoords' .* us

    #
    id1 = skipFirst ? 2 : 1
    id2 = skipLast ? 1 : 0

    node_ids = appendNodes!(mesh, nodes[id1:end-id2,:])

    return node_ids
end

function meshBoundary!(mesh::Mesh, boundary::CircularBoundary, vertices, skipFirst::Bool, skipLast::Bool)

    coords1 = vertices[boundary.vertice_start]
    coords2 = vertices[boundary.vertice_end]

    circcenter = circlecenter(coords1, coords2, boundary.radius)

    us = getFac(boundary)
    nodes = hcat(circlesection(coords1, coords2, circcenter, us)...)

    #
    id1 = skipFirst ? 2 : 1
    id2 = skipLast ? 1 : 0

    appendNodes!(mesh, nodes[id1:end-id2,:])
end

function meshBoundaries!(mesh::Mesh, meshdefinition::MeshDef)
    # create Dict for Node mapping, mapping of node_ids to boundary respectiveley vertice ids
    nodeMapping = Dict(:vertice => Dict{Int64, Int64}(), :boundary => Dict{Int64, Vector{Int64}}())

    # create nodes on bounds
    for (boundary_id, boundary) in enumerate(meshdefinition.bounds)

        # when first or last node (at vertex) is already created with other boundary -> skip
        skipFirst = skipLast = false
        if verthasNode(nodeMapping, boundary.vertice_start)
            skipFirst = true
        end

        if verthasNode(nodeMapping, boundary.vertice_end)
            skipLast = true
        end

        # create Boundary Notes and append them to mesh
        node_ids = meshBoundary!(mesh, boundary, meshdefinition.vertices, skipFirst, skipLast)

        if sign(boundary_id) == -1
            node_ids = node_ids[end:-1:1]
        end

        nodeMapping[:boundary][abs(boundary_id)] = node_ids

        # add created nodes to nodeMapping
        if !skipFirst
            nodeMapping[:vertice][boundary.vertice_start] = node_ids[1]
        else
            insert!(nodeMapping[:boundary][abs(boundary_id)], 1, nodeMapping[:vertice][boundary.vertice_start])
        end
        if !skipLast
            nodeMapping[:vertice][boundary.vertice_end] = node_ids[end]
        else
            push!(nodeMapping[:boundary][abs(boundary_id)], nodeMapping[:vertice][boundary.vertice_end])
        end
    end

    return nodeMapping
end

"""
    mesh(meshdef)

Generate mesh from meshdefinition
"""
function mesh(meshdefinition::MeshDef)

    # disable garbage collection during mesh creation
    GC.enable(false)

    # create an empty Mesh object
    mesh = emptyMesh()

    # generate nodes on boundaries
    nodeMapping = meshBoundaries!(mesh, meshdefinition)

    # generate blocks
    for (block_id, block) in enumerate(meshdefinition.blocks)
        meshBlock!(mesh, block_id, meshdefinition, nodeMapping)
    end

    # enable garbage collection again
    GC.enable(true)

    return mesh
end

function verthasNode(nodeMapping, vert_id)
    return haskey(nodeMapping[:vertice], vert_id)
end

function getBoundaryNodeIds(bound_id, nodeMapping)
    boundarynodes = nodeMapping[:boundary][abs(bound_id)]

    # sign is used as direction
    if bound_id > 0
        return boundarynodes
    else
        return boundarynodes[end:-1:1]
    end
end

function meshBlock!(mesh::Mesh, block_id, meshdef, nodeMapping)
    block = meshdef.blocks[block_id]
    if block[:type] == :transfinite
        meshTransfiniteBlock2!(mesh, block_id, meshdef, nodeMapping)
    elseif block[:type] == :transition
        meshTransitionBlock!(mesh, block, meshdef.bounds, nodeMapping)
    else
        blocktype_str = string(block[:type])
        throw(UndefVarError("Not defined blocktype ($blocktype_str)"))
    end
end

function blockdef2fun(mesh, block, bounds, nodeMapping)
    # extract blocks boundaries
    boundary_ids = block[:bounds]
    boundaries = bounds[abs.(boundary_ids)]

    # extract vertice coordinates
    firstvert(boundary_id) = boundary_id>0 ? bounds[boundary_id].vertice_start : bounds[-boundary_id].vertice_end
    vertice_ids = [firstvert(boundary_id) for boundary_id in boundary_ids]
    vertices = Dict(id=>mesh.nodes[nodeMapping[:vertice][id]] for id in vertice_ids)

    # define corner points
    P12, P14, P34, P32 = (vertices[i] for i in vertice_ids)

    # extract edge definitions
    linefunc(p1, p2) = (u) -> p1 .+ (p2.-p1) .* u
    c1 = linefunc(P12, P14)
    c2 = linefunc(P12, P32)
    c3 = linefunc(P32, P34)
    c4 = linefunc(P14, P34)

    return (P12, P14, P34, P32), (c1, c2, c3, c4)
end

"""
    blockdef2fun2(block, meshdef)

Get cornerpoints and boundary description as function
"""
function blockdef2fun2(block, meshdef)
    # extract blocks boundaries
    boundary_ids = block[:bounds]
    boundaries = meshdef.bounds[abs.(boundary_ids)]

    ps = []
    cs = []

    for i in 1:4

        p, c = boundaryfun(boundaries[i], meshdef.vertices, (sign(boundary_ids[i])==-1))

        push!(ps, p)
        push!(cs, c)
    end

    return tuple(ps...), (cs[1], u->cs[4](1-u), u->cs[3](1-u), cs[2])
end

function meshTransfiniteBlock2!(mesh::Mesh, block_id, meshdef, nodeMapping)

    elem_type = meshdef.blocks[block_id][:elem_type]

    function addElements!(mesh, last_node_ids, node_ids)
        for j in 2:size(node_ids,1)
            el_node1 = last_node_ids[j-1]
            el_node2 = node_ids[j-1]
            el_node3 = node_ids[j]
            el_node4 = last_node_ids[j]

            el_nodes = [el_node1, el_node2, el_node3, el_node4]

            el_id = appendElement!(mesh, el_nodes, elem_type)
        end
    end

    # TODO: check definition - same node_num/bias, same direction

    block = meshdef.blocks[block_id]

    # extract blocks boundaries
    boundary_ids = block[:bounds]
    boundaries = meshdef.bounds[abs.(boundary_ids)]

    # extract
    #points, functions = blockdef2fun(mesh, block, bounds, nodeMapping)
    points, functions = blockdef2fun2(block, meshdef)

    P12, P14, P34, P32 = points
    c1, c2, c3, c4 = functions

    # calculate node positions

    us = getFac(boundaries[1], flipdirection=sign(boundary_ids[1])==-1)
    vs = getFac(boundaries[2], flipdirection=sign(boundary_ids[2])==-1)

    nodearray = calculateTransfiniteNodes(c1, c2, c3, c4, P12, P14, P34, P32, us, vs)

    # add nodes and create elements
    N1 = boundaries[1].node_num
    N2 = boundaries[2].node_num

    last_node_ids = getBoundaryNodeIds(boundary_ids[4], nodeMapping)[end:-1:1]
    node_ids = zeros(Int64, size(last_node_ids, 1))

    for i in 2:N1-1
        # add bound node to node ids

        for j in 2:N2-1
            nodecoords = collect(nodearray[i, j])
            node_ids[j] = appendNodes!(mesh, [nodecoords,])[1]
        end

        # add bound node to node ids
        node_ids[1] = getBoundaryNodeIds(boundary_ids[1], nodeMapping)[i]
        node_ids[end] = getBoundaryNodeIds(boundary_ids[3], nodeMapping)[end-(i-1)]

        # add elements to mesh
        addElements!(mesh, last_node_ids, node_ids)

        last_node_ids[1:end] = node_ids[1:end]

    end

    node_ids[1:end] = getBoundaryNodeIds(boundary_ids[2], nodeMapping)

    addElements!(mesh, last_node_ids, node_ids)

end

function meshTransitionUnit!(mesh::Mesh, nodes1::Vector{Int64}, nodes2::Vector{Int64}, left::Bool, elem_type::String)
    cornerpoints_l = mesh.nodes[[nodes1[1], nodes1[3], nodes2[1], nodes2[2]]]
    cornerpoints = mesh.nodes[[nodes1[1], nodes1[end], nodes2[1], nodes2[end]]]
    cornerpoints_r = mesh.nodes[[nodes1[3], nodes1[end], nodes2[2], nodes2[end]]]

    center_ids = zeros(Int64, 3)

    for (i, cp) in enumerate((cornerpoints_l, cornerpoints, cornerpoints_r))
        centerpoint = sum(cp)/size(cp,1)

        center_ids[i] = appendNodes!(mesh, [centerpoint,])[1]
    end

    # add elements

    n1 = nodes1
    n2 = nodes2
    c = center_ids

    appendElement!(mesh, [n1[1], n1[2], c[1], n2[1]], elem_type)
    appendElement!(mesh, [n1[2], n1[3], c[2], c[1]], elem_type)
    appendElement!(mesh, [c[1], c[2], n2[2], n2[1]], elem_type)
    appendElement!(mesh, [n1[3], n1[4], c[3], c[2]], elem_type)
    appendElement!(mesh, [c[2], c[3], n2[3], n2[2]], elem_type)
    appendElement!(mesh, [n1[4], n1[5], n2[3], c[3]], elem_type)
end


text:
signatures(f) = map(m -> m.sig, methods(f).ms)
is_sf(f) = any(map(m -> Trace in m.parameters, signatures(f)))
check_sf(f) = !is_sf(f) && error("$f does not appear to be a stochastic function. (Try type annotation if it is.)")
export is_sf, check_sf

@doc raw"""
    function inference(f, method::Forward; params = (), inference_params = Dict())

No additional arguments.
"""
function inference(f, method::Forward; params = (), inference_params = Dict())
    forward_sampling(f; params = params, num_iterations = get(inference_params, "num_iterations", 1))
end

@doc raw"""
    function inference(f, method::LikelihoodWeighting; params = (), inference_params = Dict())

No additional arguments.
"""
function inference(f, method::LikelihoodWeighting; params = (), inference_params = Dict())
    likelihood_weighting(f, params...; nsamples = get(inference_params, "num_iterations", 1))
end

@doc raw"""
    inference(f; params = (), inference_params = Dict())

Do inference on the stochastic function `f`, passing as arguments to the inference algorithm
the tuple of `params` and any additional arguments to the inference algorithm in `inference_params`.
To specify a particular inference algorithm pass an `InferenceType` as the second argument
to the function, i.e.,
`inference(f, i::InferenceType; params = (), inference_params = Dict())`.
All methods accept `inference_params["num_iterations"]::Int64` as the number of iterations of
inference to perform; the interpretation of that is inference algorithm-dependent.
**This default method performs inference using `LikelihoodWeighting()`.**
"""
inference(f; params = (), inference_params = Dict()) = inference(f, LW; params = params, inference_params = inference_params)

@doc raw"""
    function inference(f, method::ImportanceSampling; params = (), inference_params = Dict())

Additional arguments: `inference_params["kernel"] = q`
where `q` is a proposal kernel for importance sampling
"""
function inference(f, method::ImportanceSampling; params = (), inference_params = Dict())
    !("kernel" in keys(inference_params)) && error("Must pass importance sampling kernel in inference_params")
    kernel = inference_params["kernel"]
    check_sf(kernel)
    importance_sampling(f, kernel; params = params, nsamples = get(inference_params, "num_iterations", 1))
end

@doc raw"""
    function inference(f, method::Nested; params = (), inference_params = Dict())

No additional arguments.
"""
function inference(f, method::Nested; params = (), inference_params = Dict())
    nested(f, rejection; params = params, num_points = get(inference_params, "num_iterations", 2))
end

@doc raw"""
    function inference(f, method::Metropolis; params = (), inference_params = Dict())

Additional arguments: `inference_params["kernels"]`
which is a vector of >= 1 proposal kernels. If this is not passed, this method uses
ancestor MH sampling.
"""
function inference(f, method::Metropolis; params = (), inference_params = Dict())
    with_kernels = "kernels" in keys(inference_params)
    if with_kernels
        kernels = inference_params["kernels"]
        map(k -> check_sf(k), kernels)
        mh(
            f, kernels;
            params = params,
            burn = get(inference_params, "burn", 0),
            thin = get(inference_params, "thin", 1),
            num_iterations = get(inference_params, "num_iterations", 2),
            inverse_verbosity = get(inference_params, "inverse_verbosity", 100)
        )
    else
        @info "No kernel(s) passed, using ancestor sampling Metropolis."
        mh(
            f;
            params = params,
            burn = get(inference_params, "burn", 0),
            thin = get(inference_params, "thin", 1),
            num_iterations = get(inference_params, "num_iterations", 2)
        )
    end
end

export inference


text:
using PaddleChainRules.Paddle: paddle, PaddleModuleWrapper, PaddleFCNet

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, DiffEqFlux
using Quadrature, Cubature, Optimisers
import ModelingToolkit: Interval, infimum, supremum

@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# 2D PDE
eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u(0,y) ~ 0.0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
       u(x,0) ~ 0.0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ Interval(0.0,1.0),
           y ∈ Interval(0.0,1.0)]

# Discretization
dx = 0.1
# number of dimensions
dim = 2

# pde system
@named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])

# callback
cb = function (p,l)
    println("Current loss is: $l")
    return false
end

# paddle
# error when setting float32
paddle.set_default_dtype("float64")
paddleNN = paddle.nn.Sequential(
        paddle.nn.Linear(2, 16),
        paddle.nn.Sigmoid(),
        paddle.nn.Linear(16, 16),
        paddle.nn.Sigmoid(),
        paddle.nn.Linear(16, 1)
    )

paddlewrap = PaddleModuleWrapper(paddleNN)

# or use a constructor for full connected network
paddlewrap = PaddleFCNet(2, 1, 3, 16;dtype="float64", activation="sigmoid")

initθ, _ = Optimisers.destructure(paddlewrap)
discretization = PhysicsInformedNN(paddlewrap, StochasticTraining(100;bcs_points = 40), init_params =initθ)
## or other training strategy
# discretization = PhysicsInformedNN(paddlewrap, GridTraining(dx), init_params =initθ)
## QuadratureTraining() works very slow here.
# discretization = PhysicsInformedNN(paddlewrap, QuadratureTraining(), init_params =initθ)
prob = discretize(pde_system,discretization)

# solve system
res = GalacticOptim.solve(prob, DiffEqFlux.ADAM(0.1); cb = cb, maxiters=4000)
prob = remake(prob, u0=res.minimizer)
res = GalacticOptim.solve(prob, DiffEqFlux.ADAM(0.01); cb = cb, maxiters=2000)
phi = discretization.phi

# error
xs,ys = [infimum(d.domain):dx/10:supremum(d.domain) for d in domains]
analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)

u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
diff_u = abs.(u_predict .- u_real)
maximum(diff_u)

max_stars_repo_name: PangHua/InfiniteScience

max_issues_repo_name: PangHua/InfiniteScience

max_forks_repo_name: PangHua/InfiniteScience


text:
#===============================================================================
  Conjugate Gradient Descent

  We implement a generic conjugate gradient descent method that includes many
  different β-choices, preconditioning, precise line searches, and more.

  A conjugate gradient method here will have to implement a CGUpdate type (CD,
  HZ, HS, ...) that controls the β update via the update_parameter function.

  The update formulae are mostly taken from HZ2006

  Let N be a neighborhood of {x ∈ R^n : f(x) ⩽ f(x₀)}

  [LA; Lipschitz Assumption]: ∃ L < ∞ : ||∇f(x)-∇f(y)|| ⩽ L||x-y|| for x,y ∈ N

  Todos:
  Might consider Daniel 1967 that uses second order
===============================================================================#

abstract type CGUpdate end
struct ConjugateGradient{Tu, TP}
  update::Tu
  P::TP
end
ConjugateGradient(;update=HZ()) = ConjugateGradient(update, nothing)
hasprecon(::ConjugateGradient{<:Any, <:Nothing}) = NoPrecon()
hasprecon(::ConjugateGradient{<:Any, <:Any}) = HasPrecon()

struct CGVars{T1, T2, T3}
    y::T1 # change in successive gradients
    d::T2 # search direction
    α::T3
    β::T3
    ls_success::Bool
end

function prepare_variables(problem, approach::LineSearch{<:ConjugateGradient, <:Any, <:Any}, x0, ∇fz)
    z = x0
    x = copy(z)
    fz, ∇fz = upto_gradient(problem, ∇fz, x)

    fx = copy(fz)
    ∇fx = copy(∇fz)

    Pg = approach.scheme.P isa Nothing ? ∇fz : copy(∇fz)
    return (x=x, fx=fx, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=nothing, Pg=Pg)
end

summary(cg::ConjugateGradient) = "Conjugate Gradient Descent ($(summary(cg.update)))"
#===============================================================================
  Conjugate Descent [Fletcher] (CD)

  R. Fletcher. "Practical Methods of Optimization vol. 1: Unconstrained
  Optimization." John Wiley & Sons, New York (1987).
===============================================================================#
struct CD <: CGUpdate end
function update_parameter(mstyle, cg::CD, d, ∇fz, ∇fx, y, P, P∇fz)
  -dot(∇fz, P∇fz)/dot(d, ∇fx)
end

#===============================================================================
  Hager-Zhang (HZ)
===============================================================================#
struct HZ{Tη} <: CGUpdate
    η::Tη # a "forcing term"
end
HZ() = HZ(0.4)
function update_parameter(mstyle, cg::HZ, d, ∇fz, ∇fx, y, P, P∇fz)
    T = eltype(∇fz)
    θ = T(2)
    η = T(cg.η)
    # βHS = update_parameter(mstyle, HS(), d,  ∇fz, ∇fx, y, P, P∇fz)
    # but let's save the dy calculation from begin repeated for
    # efficiency's sake
    dy = dot(d, y)
    βHS = dot(y, P∇fz)/dy

    # Apply preconditioner to y
    Py = apply_preconditioner(mstyle, P, copy(P∇fz), y)
    βN = βHS - θ*dot(y, Py)/dy*dot(d, ∇fz)/dy

    # 2013 version is scale invariant
    Pinv_y = apply_inverse_preconditioner(mstyle, P, copy(P∇fz), y)
    ηk = η*dot(d, ∇fx)/dot(d, Pinv_y)
    # 2006 version
    # ηk = -inv(norm(d)*min(T(cg.η), norm(∇fx)))

    βkp1 = max(βN, ηk)
end

#===============================================================================
  Hestenes-Stiefel (HS)

  Jamming protection from y in the numerator. x-x' small means ∇f(x)-∇f(x')
  should also be small.

  Hestenes, Magnus R., and Eduard Stiefel. "Methods of conjugate gradients for
  solving linear systems." Journal of research of the National Bureau of
  Standards 49, no. 6 (1952): 409-436.
===============================================================================#
struct HS <: CGUpdate end
function update_parameter(mstyle, cg::HS, d, ∇fz, ∇fx, y, P, P∇fz)
  dot(y, P∇fz)/dot(d, y)
end
#===============================================================================
  Fletcher-Reeves (FR)

  Fletcher, Reeves, and Colin M. Reeves. "Function minimization by conjugate
  gradients." The computer journal 7, no. 2 (1964): 149-154.
===============================================================================#
struct FR <: CGUpdate end
function update_parameter(mstyle, cg::FR, d, ∇fz, ∇fx, y, P, P∇fz)
    dot(∇fz, P∇fz)/dot(∇fx, ∇fx)
end
#===============================================================================
  Polak-Ribiére-Polyak (PRP)

  Jamming protection from y in the numerator. x-x' small means ∇f(x)-∇f(x')
  should also be small.

  Polak, Elijah, and Gerard Ribiere. "Note sur la convergence de méthodes de
  directions conjuguées." ESAIM: Mathematical Modelling and Numerical Analysis-
  Modélisation Mathématique et Analyse Numérique 3, no. R1 (1969): 35-43.

  Polyak, Boris T. "The conjugate gradient method in extremal problems." USSR
  Computational Mathematics and Mathematical Physics 9, no. 4 (1969): 94-112.
===============================================================================#
struct PRP{Plus} end
PRP(;plus=true) = PRP{plus}()
function update_parameter(mstyle, ::PRP, d, ∇fz, ∇fx, y, P, P∇fz)
    dot(y, P∇fz)/dot(∇fx, ∇fx)
end
function update_parameter(::PRP{true}, d, ∇fz, ∇fx, y, P, P∇fz)
  βPR = update_parameter(RP{false}(), d, ∇fz, ∇fx, y, P, P∇fz)
  max(0, β)
end

#===============================================================================
  Liu-Storey (LS)

  Liu-Storey is identical to PRP for exact line search.

  Jamming protection from y in the numerator. x-x' small means ∇f(x)-∇f(x')
  should also be small.

  Liu, Y., and C. Storey. "Efficient generalized conjugate gradient algorithms,
  part 1: theory." Journal of optimization theory and applications 69, no. 1
  (1991): 129-137.
===============================================================================#
struct LS <: CGUpdate end
function update_parameter(mstyle, ::LS, d, ∇fz, ∇fx, y, P, P∇fz)
  -dot(y, P∇fz)/dot(d, ∇fx)
end
#===============================================================================
  Dai-Yuan (DY)

  Different from CD and FR. With Wolfe line search it always generates descent
  directions. Global convergence with Lipchitz Assumption [LA].

  Dai, Yu-Hong, and Yaxiang Yuan. "A nonlinear conjugate gradient method with a
  strong global convergence property." SIAM Journal on optimization 10, no. 1
  (1999): 177-182.
===============================================================================#
struct DY <: CGUpdate end
function update_parameter(mstyle, cg::DY, d, ∇fz, ∇fx, y, P, P∇fz)
  dot(∇fz, P∇fz)/dot(d, y)
end

#===============================================================================
  Wei-Yao-Liu (VPRP)

  Yu, Gaohang, Yanlin Zhao, and Zengxin Wei. "A descent nonlinear conjugate
  gradient method for large-scale unconstrained optimization." Applied
  mathematics and computation 187, no. 2 (2007): 636-643.
===============================================================================#
struct VPRP <: CGUpdate end
function update_parameter(mstyle, cg::VPRP, d, ∇fz, ∇fx, y, P, P∇fz)
  a = dot(∇fz, P∇fz)
  b = dot(∇fx, P∇fz)
  c = norm(∇fx, 2)
  num = a - sqrt(a)/c*b
  num/dot(d, y)
end

# using an "initial vectors" function we can initialize s if necessary or nothing if not to save on vectors
function solve(problem::OptimizationProblem, x0, cg::ConjugateGradient, options::OptimizationOptions)
  _solve(problem, x0, LineSearch(cg, HZAW()), options, mstyle(problem))
end
function solve(problem::OptimizationProblem, x0, approach::LineSearch{<:ConjugateGradient, <:LineSearcher}, options::OptimizationOptions)
  _solve(problem, x0, approach, options, mstyle(problem))
end
function _solve(problem::OptimizationProblem, x0, approach::LineSearch{<:ConjugateGradient, <:LineSearcher}, options::OptimizationOptions, mstyle::MutateStyle)
    t0 = time()
    #==============
         Setup
    ==============#
    Tx = eltype(x0)

    objvars = prepare_variables(problem, approach, x0, copy(x0))
    f0, ∇f0 = objvars.fx, norm(objvars.∇fx, Inf) # use user norm

    y, d, α, β = copy(objvars.∇fz), -copy(objvars.∇fx), Tx(0), Tx(0)
    cgvars = CGVars(y, d, α, β, true)

    k = 1
    objvars, P, cgvars = iterate(mstyle, cgvars, objvars, approach, problem, options)
    is_converged = converged(approach, objvars, ∇f0, options)
    while k < options.maxiter && !any(is_converged)
        k += 1
        objvars, P, cgvars = iterate(mstyle, cgvars, objvars, approach, problem, options, P, false)
        is_converged = converged(approach, objvars, ∇f0, options)
    end
    x, fx, ∇fx, z, fz, ∇fz, B = objvars
    return ConvergenceInfo(approach, (beta=β, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fx=∇fx, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=k, time=time()-t0), options)
end
function iterate(mstyle::InPlace, cgvars::CGVars, objvars, approach::LineSearch{<:ConjugateGradient, <:Any, <:Any}, problem::OptimizationProblem, options::OptimizationOptions, P=nothing, is_first=nothing)
    # split up the approach into the hessian approximation scheme and line search
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    Tx = eltype(x)

    scheme, linesearch = modelscheme(approach), algorithm(approach)
    y, d, α, β = cgvars.y, cgvars.d, cgvars.α, cgvars.β

    # Move nexts into currs
    fx = fz
    copyto!(x, z)
    copyto!(∇fx, ∇fz)

    # Precondition current gradient and calculate the search direction
    P = update_preconditioner(scheme, x, P)
    P∇fz = apply_preconditioner(mstyle, P, Pg, ∇fz)
    @. d = -P∇fz + β*d

    α_0 = initial(scheme.update, a->value(problem, x.+a.*d), α, x, fx, dot(d, ∇fx), ∇fx, is_first)
    φ = _lineobjective(mstyle, problem, ∇fz, z, x, d, fx, dot(∇fx, d))

    # Perform line search along d
    α, f_α, ls_success = find_steplength(mstyle, linesearch, φ, Tx(1))

    # Calculate final step vector and update the state
    if ls_success
        z = retract(problem, z, x, d, α)
        fz, ∇fz = upto_gradient(problem, ∇fz, z)
        @. y = ∇fz - ∇fx
    else
        # if no succesful search direction is found, reset to gradient
        y .= .-∇fz
    end
    β = update_parameter(mstyle, scheme.update, d, ∇fz, ∇fx, y, P, P∇fz)

    return (x=x, fx=fx, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=nothing, Pg=Pg), P, CGVars(y, d, α, β, ls_success)
end

function iterate(mstyle::OutOfPlace, cgvars::CGVars, objvars, approach::LineSearch{<:ConjugateGradient, <:Any, <:Any}, problem::OptimizationProblem, options::OptimizationOptions, P=nothing, is_first=nothing)
    # split up the approach into the hessian approximation scheme and line search
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    Tx = eltype(x)
    scheme, linesearch = modelscheme(approach), algorithm(approach)
    d, α, β = cgvars.d, cgvars.α, cgvars.β

    # Move nexts into currs
    fx = fz
    x = copy(z)
    ∇fx = copy(∇fz)

    # Update preconditioner
    P = update_preconditioner(scheme, x, P)
    P∇fz = apply_preconditioner(mstyle, P, Pg, ∇fz)
    # Update current gradient and calculate the search direction
    d = @. -P∇fz + β*d

    α_0 = initial(scheme.update, a-> value(problem, x .+ a.*d), α, x, fx, dot(d, ∇fx), ∇fx, is_first)
    φ = _lineobjective(mstyle, problem, ∇fz, z, x, d, fx, dot(∇fx, d))

    # Perform line search along d
    α, f_α, ls_success = find_steplength(mstyle, linesearch, φ, Tx(1))

    z = retract(problem, z, x, d, α)
    fz, ∇fz = upto_gradient(problem, ∇fz, z)
    y = @. ∇fz - ∇fx
    β = update_parameter(mstyle, scheme.update, d, ∇fz, ∇fx, y, P, P∇fz)
    return (x=x, fx=fx, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=nothing, Pg=Pg), P, CGVars(y, d, α, β, ls_success)
end

function initial(cg, φ, α, x, φ₀, dφ₀, ∇fx, is_first)
    T = eltype(x)
    ψ₀ = T(0.01)
    ψ₁ = T(0.1)
    ψ₂ = T(2.0)
    quadstep = true
    if is_first isa Nothing
        if !all(x .≈ T(0)) # should we define "how approx we want?"
            return ψ₀ * norm(x, Inf)/norm(∇fx, Inf)
        elseif !(φ₀ ≈ T(0))
            return ψ₀ * abs(φ₀)/norm(∇fx, 2)^2
        else
            return T(1)
        end
    elseif quadstep
        R = ψ₁*α
        φR = φ(R)
        if φR ≤ φ₀
            c = (φR - φ₀ - dφ₀*R)/R^2
            if c > 0
               return -dφ₀/(T(2)*c) # > 0 by df0 < 0 and c > 0
            end
        end
    end
    return ψ₂*α
end


text:
function irf_ci_asymptotic(V::VAR, H)
    (T,K) = size(V.Y')
    if V.i == true
        SIGa = kron(inv(V.X*V.X'/(T-V.p)),V.Σ)
        SIGa = SIGa[K+1:end,K+1:end]
    else
        SIGa = kron(inv(V.X*V.X'/(T-V.p)),V.Σ)
    end
    # Calculation of stdev follows Lutkepohl(2005) p.111,93
    A = get_VAR1_rep(V)
    A0inv = full(chol(V.Σ)')
    STD   = zeros(K^2,H+1)
    COV2   = zeros(K^2,H+1)
    J = [eye(K) zeros(K,K*(V.p-1))]
    L = elimat(K)
    Kk = commutation(K,K)
    Hk = L'/(L*(eye(K^2)+Kk)*kron(A0inv,eye(K))*L')
    D = duplication(K)
    Dplus = (D'*D)\D'
    SIGsig = 2*Dplus*kron(V.Σ,V.Σ)*Dplus';
    Cbar0 = kron(eye(K),J*eye(K*V.p)*J')*Hk;
    STD[:,1] = vec((reshape(diag(real(sqrt(complex(Cbar0*SIGsig*Cbar0'/(T-V.p))))),K,K))')
    COV2[:,1] = vec((reshape(diag((Cbar0*SIGsig*Cbar0'/(T-V.p))),K,K))');
    for h=1:H
        Gi = zeros(K^2,K^2*V.p)
        for m=0:(h-1)
            Gi += kron(J*(A')^(h-1-m),J*(A^m)*J')
        end
        C = kron(A0inv',eye(K))*Gi
        Cbar = kron(eye(K),J*A^h*J')*Hk
        STD[:,h+1] = vec((reshape(diag(real(sqrt(complex(C*SIGa*C'+Cbar*SIGsig*Cbar')/(T-V.p)))),K,K))')
        COV2[:,h+1] = vec((reshape(diag(((Cbar*SIGsig*Cbar')/(T-V.p))),K,K))')
    end
    return STD,COV2
end


text:
# specialized likelihood for MCMC using net
function LL(θ, m, S, model, info)
    k = size(m,1)
    ms = zeros(S, k)
    Threads.@threads for s = 1:S
        ms[s,:] = Float64.(model(transform(WileE_model(θ[:])', info)'))
    end
    mbar = mean(ms,dims=1)[:]
    Σ = cov(ms)
    x = (m .- mbar)[:]
    lnL = try
        lnL =-0.5*log(det(Σ)) - 0.5*x'*inv(Σ)*x # for Bayesian
    catch
        lnL = -Inf
    end
    return lnL
end

# version without net
function LL(θ, m, S)
    k = size(m,1)
    ms = zeros(S, k)
    Threads.@threads for s = 1:S
        ms[s,:] = WileE_model(θ[:])
    end
    mbar = mean(ms,dims=1)[:]
    Σ = cov(ms)
    x = (m .- mbar)[:]
    lnL = try
        lnL =-0.5*log(det(Σ)) - 0.5*x'*inv(Σ)*x # for Bayesian
    catch
        lnL = -Inf
    end
    return lnL
end


text:
"""
Program to check if a number is an [Armstrong/Narcissistic number](https://en.wikipedia.org/wiki/Narcissistic_number) in decimal system.

Armstrong number is a number that is the sum of its own digits raised to the power of the number of digits.

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)

A positive integer is called an Armstrong number (of order n) if

`abcd... = a^n + b^n + c^n + d^n +....`
"""
function is_armstrong(x)
    # If the number is less than 0, then it is not a armstrong number.
    if x < 0
        return false
    elseif x == 0
        return true
    else
        # Finding the total number of digits in the number
        length = Int(floor(log10(x) + 1))
        result = 0
        temp = x
        for i in 1:length
            val = temp % 10
            # Finding each digit raised to the power total digit and add it to the total sum
            result = result + Int(trunc(val))^length
            temp = temp / 10
        end
        return x == result ? true : false
    end
end

max_stars_repo_name: madisonestabrook/Julia

max_issues_repo_name: madisonestabrook/Julia

max_forks_repo_name: madisonestabrook/Julia


text:
struct QueryableSource <: Queryable
    getiterator
end


text:
using GVFHordes
import ProgressMeter

StatsBase.sample(rng::Random.AbstractRNG, p::GVFHordes.GVFParamFuncs.FunctionalPolicy, s, actions) =
    sample(rng, Weights([p.func(state_t = s, action_t = a) for a in actions]))

"""
    monte_carlo_return

Perform a monte_carlo_return using the provided environment, gvf, and start_state. For an environment to be compatible it must
implement `reset!(env, state)`

The interface is limited in that it can't do rollouts on compositional GVFs.
"""
function monte_carlo_return(env,
                            gvf,
                            start_state,
                            start_action,
                            num_returns,
                            γ_thresh=1e-6,
                            max_steps=Int(1e6),
                            rng=Random.GLOBAL_RNG)
    returns = zeros(num_returns)

    for ret in 1:num_returns
        step = 0
        term = false
        cumulative_gamma = 1.0

        cur_state = start!(env, start_state)
        # next_action = StatsBase.sample(rng, GVFHordes.policy(gvf), cur_state, get_actions(env))
        next_action = start_action

        while cumulative_gamma > γ_thresh &&
            step < max_steps &&
            term == false
            # Take action
            action = next_action
            next_state, r, term = MinimalRLCore.step!(env, action, rng)

            # Get next action for GVFs
            next_action = StatsBase.sample(rng, GVFHordes.policy(gvf), next_state, get_actions(env))

            # Update Return
            # c, γ, pi_prob = get(gvf, cur_state, action, next_state, next_action, nothing)
            c, γ, pi_prob = get(gvf; state_t = cur_state, action_t = action, state_tp1 = next_state, action_tp1 = next_action)
            returns[ret] += cumulative_gamma*c
            cumulative_gamma *= γ*(1-term)

            cur_state = next_state
            step += 1
        end
    end

    return returns
end

function monte_carlo_returns(env, gvf, start_states, actions, num_returns, γ_thresh, max_steps=Int(1e6), rng=Random.GLOBAL_RNG)

    states_actions = zip(start_states, actions)
    ret = Vector{Vector{Float64}}(undef, length(states_actions))
    prg_meter = ProgressMeter.Progress(length(states_actions))
    Threads.@threads for i ∈ 1:length(states_actions)
        t_env = deepcopy(env)
        ret[i] = monte_carlo_return(t_env, gvf, start_states[i], actions[i], num_returns, γ_thresh, max_steps)
        ProgressMeter.next!(prg_meter)
    end
    ret

end

function monte_carlo_returns_agg(env, gvf, start_states, actions, num_returns, γ_thresh, max_steps=Int(1e6), rng=Random.GLOBAL_RNG; agg)

    # states_actions = zip(start_states, actions)
    ret = Vector{Float64}(undef, length(start_states))
    prg_meter = ProgressMeter.Progress(length(start_states))

    Threads.@threads for i ∈ 1:length(start_states)
        t_env = deepcopy(env)
        ret[i] = agg(monte_carlo_return(t_env, gvf, start_states[i], actions[i], num_returns, γ_thresh, max_steps))
        ProgressMeter.next!(prg_meter)
    end
    ret

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using GLAbstraction, GLPlot, Reactive
GLPlot.init()
text = "whatup internet!?\n#Crusont"

glplot(text, color = [RGBA(rand(), rand(), 0,1) for i=1:length(text)]) # You can either, supply a texture with colors


text:
using Plots
x = 1:10
y = rand(10)
plot(x, y)


text:
import Literate

function preprocess_example_shared(
        content::AbstractString;
        output_directory::AbstractString,
        )::String
    content = replace(content,
                      "## %PREDICTMD_GENERATED_BY%\n" =>
                          string("## This file was generated by PredictMD ",
                                 "version $(version())",
                                 "\n",
                                 "## For help, please visit ",
                                 "https://predictmd.net",
                                 "\n",))
    return content
end

function get_if_include_test_statements_regex()::Regex
    test_statements_regex::Regex = r"# PREDICTMD IF INCLUDE TEST STATEMENTS\n([\S\s]*?)# PREDICTMD ELSE\n([\S\s]*?)# PREDICTMD ENDIF INCLUDE TEST STATEMENTS\n{0,5}"
    return test_statements_regex
end

function preprocess_example_do_not_include_test_statements(
        content::AbstractString;
        output_directory::AbstractString,
        )::String
    content::String = preprocess_example_shared(
        content;
        output_directory = output_directory,
        )
    test_statements_regex::Regex = get_if_include_test_statements_regex()
    for m in eachmatch(test_statements_regex, content)
        original_text::String = String(m.match)
        replacement_text = string(strip(String(m[2])), "\n\n",)
        content = replace(content, original_text => replacement_text)
    end
    content = replace(content, r"logger_stream, " => "")
    return content
end

function preprocess_example_include_test_statements(
        content::AbstractString;
        output_directory::AbstractString,
        )::String
    content::String = preprocess_example_shared(
        content;
        output_directory = output_directory,
        )
    test_statements_regex::Regex = get_if_include_test_statements_regex()
    for m in eachmatch(test_statements_regex, content)
        content = replace(content,
                          String(m.match) =>
                              string(strip(String(m[1])), "\n\n",))
    end
    return content
end

function fix_example_blocks(filename::AbstractString)::Nothing
    if filename_extension(filename) == ".md"
        content = read(filename, String)
        rm(filename; force = true, recursive = true,)
        content = replace(content, r"```@example \w*\n" => "```julia\n")
        write(filename, content)
    end
    return nothing
end

function generate_examples(
        output_directory::AbstractString;
        execute_notebooks = false,
        markdown = false,
        notebooks = false,
        scripts = false,
        include_test_statements::Bool = false,
        )::String
    _abspath_output_directory::String = abspath(output_directory)
    if Sys.iswindows()
        execute_notebooks = false
    end
    ENV["PREDICTMD_IS_MAKE_EXAMPLES"] = "true"
    if !markdown && !notebooks && !scripts
        error(
            string(
                "At least one of markdown, notebooks, scripts must be true.",
                )
            )
    end
    if ispath(output_directory)
        error(
            string(
                "The output directory already exists. ",
                "Delete the output directory and then ",
                "re-run generate_examples."
                )
            )
    end

    if include_test_statements
        preprocess_example = (x) ->
            preprocess_example_include_test_statements(
                x;
                output_directory = _abspath_output_directory,
                )
    else
        preprocess_example = (x) ->
            preprocess_example_do_not_include_test_statements(
                x;
                output_directory = _abspath_output_directory,
                )
    end

    @debug("Starting to generate examples...")

    temp_examples_dir = joinpath(
        maketempdir(),
        "generate_examples",
        "PredictMDTemp",
        "docs",
        "src",
        "examples",
        )
    try
        mkpath(temp_examples_dir)
    catch
    end

    examples_input_parent_directory =
        PredictMD.package_directory("templates", "examples")

    cpu_examples_input_parent_directory = joinpath(
        examples_input_parent_directory,
        "cpu_examples",
        )
    cpu_examples_output_parent_directory = joinpath(
        temp_examples_dir,
        "cpu_examples",
        )
    try
        mkpath(cpu_examples_output_parent_directory)
    catch
    end

    boston_housing_input_directory = joinpath(
        cpu_examples_input_parent_directory,
        "boston_housing",
        )
    boston_housing_output_directory = joinpath(
        cpu_examples_output_parent_directory,
        "boston_housing",
        )

    boston_housing_input_src_directory = joinpath(
        boston_housing_input_directory,
        "src",
        )
    boston_housing_output_src_directory = joinpath(
        boston_housing_output_directory,
        "src",
        )
    try
        mkpath(boston_housing_output_directory)
        mkpath(boston_housing_output_src_directory)
    catch
    end
    for x in [".gitignore", "Project.toml", "README.md"]
        cp(joinpath(boston_housing_input_directory, x),
           joinpath(boston_housing_output_directory, x);
           force=true)
    end
    boston_housing_input_file_list =
        readdir(boston_housing_input_src_directory)
    boston_housing_input_file_list =
        boston_housing_input_file_list[
            [endswith(x, ".jl") for x in
                boston_housing_input_file_list]
            ]
    sort!(boston_housing_input_file_list)
    for input_file in boston_housing_input_file_list
        input_file_full_path = joinpath(
            boston_housing_input_src_directory,
            input_file,
            )
        if markdown
            Literate.markdown(
                input_file_full_path,
                boston_housing_output_src_directory;
                codefence = "```@example boston_housing" => "```",
                documenter = true,
                preprocess = preprocess_example,
                )
        end
        if notebooks
            Literate.notebook(
                input_file_full_path,
                boston_housing_output_src_directory;
                documenter = true,
                execute = execute_notebooks,
                preprocess = preprocess_example,
                )
        end
        if scripts
            Literate.script(
                input_file_full_path,
                boston_housing_output_src_directory;
                documenter = true,
                keep_comments = true,
                preprocess = preprocess_example,
                )
        end
    end

    breast_cancer_biopsy_input_directory = joinpath(
        cpu_examples_input_parent_directory,
        "breast_cancer_biopsy",
        )
    breast_cancer_biopsy_output_directory = joinpath(
        cpu_examples_output_parent_directory,
        "breast_cancer_biopsy",
        )

    breast_cancer_biopsy_input_src_directory = joinpath(
        breast_cancer_biopsy_input_directory,
        "src",
        )
    breast_cancer_biopsy_output_src_directory = joinpath(
        breast_cancer_biopsy_output_directory,
        "src",
        )
    try
        mkpath(breast_cancer_biopsy_output_directory)
        mkpath(breast_cancer_biopsy_output_src_directory)
    catch
    end
    for x in [".gitignore", "Project.toml", "README.md"]
        cp(joinpath(breast_cancer_biopsy_input_directory, x),
           joinpath(breast_cancer_biopsy_output_directory, x);
           force=true)
    end
    breast_cancer_biopsy_input_file_list =
        readdir(breast_cancer_biopsy_input_src_directory)
    breast_cancer_biopsy_input_file_list =
        breast_cancer_biopsy_input_file_list[
            [endswith(x, ".jl") for x in
                breast_cancer_biopsy_input_file_list]
            ]
    sort!(breast_cancer_biopsy_input_file_list)
    for input_file in breast_cancer_biopsy_input_file_list
        input_file_full_path = joinpath(
            breast_cancer_biopsy_input_src_directory,
            input_file,
            )
        if markdown
            Literate.markdown(
                input_file_full_path,
                breast_cancer_biopsy_output_src_directory;
                codefence = "```@example breast_cancer_biopsy" => "```",
                documenter = true,
                preprocess = preprocess_example,
                )
        end
        if notebooks
            Literate.notebook(
                input_file_full_path,
                breast_cancer_biopsy_output_src_directory;
                documenter = true,
                execute = execute_notebooks,
                preprocess = preprocess_example,
                )
        end
        if scripts
            Literate.script(
                input_file_full_path,
                breast_cancer_biopsy_output_src_directory;
                documenter = true,
                keep_comments = true,
                preprocess = preprocess_example,
                )
        end
    end

    for (root, dirs, files) in walkdir(temp_examples_dir)
        for f in files
            filename = joinpath(root, f)
            fix_example_blocks(filename)
        end
    end

    try
        mkpath(dirname(output_directory))
    catch
    end

    cp(
        temp_examples_dir,
        output_directory;
        force = true,
        )

    @debug(
        string(
            "Finished generating examples. ",
            "Files were written to: \"",
            output_directory,
            "\".",
            )
        )
    ENV["PREDICTMD_IS_MAKE_EXAMPLES"] = "false"
    return output_directory
end


text:
using Test, MadNLP, MadNLPHSL, MadNLPTests

testset = [
    [
        "HSL-Ma27",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPMa27,
            print_level=MadNLP.ERROR),
        []
    ],
    [
        "HSL-Ma57",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPMa57,
            print_level=MadNLP.ERROR),
        []
    ],
    [
        "HSL-Ma77",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPMa77,
            print_level=MadNLP.ERROR),
        ["unbounded"]
    ],
    [
        "HSL-Ma86",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPMa86,
            print_level=MadNLP.ERROR),
        []
    ],
    [
        "HSL-Ma97",
        ()->MadNLP.Optimizer(
            linear_solver=MadNLPMa97,
            print_level=MadNLP.ERROR),
        []
    ],
]

@testset "MadNLPHSL test" begin
    for (name,optimizer_constructor,exclude) in testset
        test_madnlp(name,optimizer_constructor,exclude)
    end
end


text:
include("utils.jl")

@testset "Euclidean" begin
    E = Euclidean(3)
    Ec = Euclidean(3; field = ℂ)
    EM = Manifolds.MetricManifold(E, Manifolds.EuclideanMetric())
    @test repr(E) == "Euclidean(3; field = ℝ)"
    @test repr(Ec) == "Euclidean(3; field = ℂ)"
    @test repr(Euclidean(2, 3; field = ℍ)) == "Euclidean(2, 3; field = ℍ)"
    @test Manifolds.allocation_promotion_function(Ec, get_vector, ()) === complex
    @test is_default_metric(EM)
    @test is_default_metric(E, Manifolds.EuclideanMetric())
    @test Manifolds.default_metric_dispatch(E, Manifolds.EuclideanMetric()) === Val{true}()
    p = zeros(3)
    @test det_local_metric(EM, p) == one(eltype(p))
    @test log_local_metric_density(EM, p) == zero(eltype(p))
    @test project!(E, p, p) == p
    @test embed!(E, p, p) == p
    @test manifold_dimension(Ec) == 2 * manifold_dimension(E)
    X = zeros(3)
    X[1] = 1.0
    Y = similar(X)
    project!(E, Y, p, X)
    @test Y == X

    # real manifold does not allow complex values
    @test_throws DomainError is_manifold_point(Ec, [:a, :b, :b], true)
    @test_throws DomainError is_manifold_point(E, [1.0, 1.0im, 0.0], true)
    @test_throws DomainError is_manifold_point(E, [1], true)
    @test_throws DomainError is_tangent_vector(Ec, [:a, :b, :b], [1.0, 1.0, 0.0], true)
    @test_throws DomainError is_tangent_vector(E, [1.0, 1.0im, 0.0], [1.0, 1.0, 0.0], true) # real manifold does not allow complex values
    @test_throws DomainError is_tangent_vector(E, [1], [1.0, 1.0, 0.0], true)
    @test_throws DomainError is_tangent_vector(E, [0.0, 0.0, 0.0], [1.0], true)
    @test_throws DomainError is_tangent_vector(E, [0.0, 0.0, 0.0], [1.0, 0.0, 1.0im], true)
    @test_throws DomainError is_tangent_vector(Ec, [0.0, 0.0, 0.0], [:a, :b, :c], true)

    @test E^2 === Euclidean(3, 2)
    @test ^(E, 2) === Euclidean(3, 2)
    @test E^(2,) === Euclidean(3, 2)
    @test Ec^(4, 5) === Euclidean(3, 4, 5; field = ℂ)

    manifolds = [E, EM, Ec]
    types = [Vector{Float64}]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    TEST_DOUBLE64 && push!(types, Vector{Double64})
    TEST_STATIC_SIZED && push!(types, MVector{3,Float64})

    types_complex = [Vector{ComplexF64}]
    TEST_FLOAT32 && push!(types_complex, Vector{ComplexF32})
    TEST_DOUBLE64 && push!(types_complex, Vector{ComplexDF64})
    TEST_STATIC_SIZED && push!(types_complex, MVector{3,ComplexF64})

    for M in manifolds
        basis_types = if M == E
            (
                DefaultOrthonormalBasis(),
                ProjectedOrthonormalBasis(:svd),
                DiagonalizingOrthonormalBasis([1.0, 2.0, 3.0]),
            )
        elseif M == Ec
            (
                DefaultOrthonormalBasis(),
                DefaultOrthonormalBasis(ℂ),
                DiagonalizingOrthonormalBasis([1.0, 2.0, 3.0]),
            )
        else
            ()
        end
        for T in types
            @testset "$M Type $T" begin
                pts = [
                    convert(T, [1.0, 0.0, 0.0]),
                    convert(T, [0.0, 1.0, 0.0]),
                    convert(T, [0.0, 0.0, 1.0]),
                ]
                test_manifold(
                    M,
                    pts,
                    test_reverse_diff = isa(T, Vector),
                    test_project_point = true,
                    test_project_tangent = true,
                    test_musical_isomorphisms = true,
                    test_default_vector_transport = true,
                    vector_transport_methods = [
                        ParallelTransport(),
                        SchildsLadderTransport(),
                        PoleLadderTransport(),
                    ],
                    test_mutating_rand = isa(T, Vector),
                    point_distributions = [Manifolds.projected_distribution(
                        M,
                        Distributions.MvNormal(zero(pts[1]), 1.0),
                    )],
                    tvector_distributions = [Manifolds.normal_tvector_distribution(
                        M,
                        pts[1],
                        1.0,
                    )],
                    basis_types_vecs = basis_types,
                    basis_types_to_from = basis_types,
                    basis_has_specialized_diagonalizing_get = true,
                    test_vee_hat = isa(M, Euclidean),
                )
            end
        end
    end
    for T in types_complex
        @testset "Complex Euclidean, type $T" begin
            pts = [
                convert(T, [1.0im, -1.0im, 1.0]),
                convert(T, [0.0, 1.0, 1.0im]),
                convert(T, [0.0, 0.0, 1.0]),
            ]
            test_manifold(
                Ec,
                pts,
                test_reverse_diff = isa(T, Vector),
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_default_vector_transport = true,
                test_vee_hat = false,
            )
        end
    end

    @testset "hat/vee" begin
        E = Euclidean(3, 2)
        p = collect(reshape(1.0:6.0, (3, 2)))
        X = collect(reshape(7.0:12.0, (3, 2)))
        @test hat(E, p, vec(X)) ≈ X
        Y = allocate(X)
        @test hat!(E, Y, p, vec(X)) === Y
        @test Y ≈ X
        @test vee(E, p, X) ≈ vec(X)
        Y = allocate(vec(X))
        @test vee!(E, Y, p, X) === Y
        @test Y ≈ vec(X)
    end

    @testset "Number systems power" begin
        @test ℝ^2 === Euclidean(2)
        @test ℝ^(2, 3) === Euclidean(2, 3)

        @test ℂ^2 === Euclidean(2; field = ℂ)
        @test ℂ^(2, 3) === Euclidean(2, 3; field = ℂ)

        @test ℍ^2 === Euclidean(2; field = ℍ)
        @test ℍ^(2, 3) === Euclidean(2, 3; field = ℍ)
    end

    @testset "Embeddings into larger Euclidean Manifolds" begin
        M = Euclidean(3, 3)
        N = Euclidean(4, 4)
        O = EmbeddedManifold(M, N)
        # first test with same length of sizes
        p = ones(3, 3)
        q = zeros(4, 4)
        qT = zeros(4, 4)
        qT[1:3, 1:3] .= 1.0
        embed!(O, q, p)
        @test norm(qT - q) == 0
        qM = embed(O, p)
        @test norm(project(O, qM) - p) == 0
        @test norm(qT - qM) == 0
        # test with different sizes, check that it only fills first element
        q2 = zeros(4, 4, 3)
        q2T = zeros(4, 4, 3)
        q2T[1:3, 1:3, 1] .= 1.0
        embed!(O, q2, p)
        @test norm(q2T - q2) == 0
        O2 = EmbeddedManifold(M, Euclidean(4, 4, 3))
        q2M = embed(O2, p)
        @test norm(q2T - q2M) == 0
        # wrong size error checks
        @test_throws DomainError embed!(O, zeros(3, 3), zeros(3, 3, 5))
        @test_throws DomainError embed!(O, zeros(3, 3), zeros(4, 4))
        @test_throws DomainError project!(O, zeros(3, 3, 5), zeros(3, 3))
        @test_throws DomainError project!(O, zeros(4, 4), zeros(3, 3))
    end

    @testset "Embedding Real into Complex" begin
        M = Euclidean(3, 3)
        N = Euclidean(3, 4; field = ℂ)
        O = EmbeddedManifold(M, N)
        p = ones(3, 3)
        qT = zeros(ComplexF64, 3, 4)
        qT[1:3, 1:3] .= 1.0
        q = embed(O, p)
        @test is_manifold_point(N, q)
        @test q == qT
        q2 = zeros(ComplexF64, 3, 4)
        embed!(O, q2, p)
        @test q2 == qT
    end
end