text:
abstract type AbstractThermoModel end #base model
abstract type AbstractHelmholtzModel <:  AbstractThermoModel end
abstract type AbstractMixingRule end # for mixing rules
abstract type AbstractPhase end #for solvers

const R_GAS = ustrip(Unitful.R)
Base.broadcastable(x::AbstractHelmholtzModel) = Ref(x)

"""
core_helmholtz(model,V,T,x)

calculates the helmholtz energy. uses the following units:
v:molar volume, m3/mol ,SI
T:Temperature, K , SI
x:molar fraction, adimensional

"""

function core_helmholtz(model::M,V,T,x) where M <:AbstractHelmholtzModel
end
# partial_molar_property
#is equivalent to df/dni with with T,v, nj constant

function partial_molar_property(core_fn,model,v,T,x)
    function partial_fun(n)
        sum_n = sum(n)
        xx = n  .* inv(sum_n)
        return sum_n*core_fn(model,v/sum_n,T,xx)
    end
    return ForwardDiff.gradient(partial_fun,x)
end

function core_logfugacity_coefficient(model,v,T,x)
    dAdn =  partial_molar_property(core_residual_helmholtz,model,v,T,x)
    lnz = log(core_compressibility_factor(model,v,T,x))
    R = R_GAS
    return dAdn/(R*T) .- lnz
end

function core_logfugacity_coefficient(model,P,v,T,x)
    dAdn =  partial_molar_property(core_residual_helmholtz,model,v,T,x)
    lnz = log(core_compressibility_factor(model,P,v,T,x))
    R = R_GAS
    return dAdn/(R*T) .- lnz
end

function core_logfugacity_coefficientn(model,v,T,n)
    return core_fugacity_coefficient(model,v,T,n/sum(n))
end

function core_logfugacity_coefficientn(model,P,v,T,n)
    return core_fugacity_coefficient(model,P,v,T,n/sum(n))
end

function core_logfugacity(model,v,T,x)
    dAdn =  partial_molar_property(core_residual_helmholtz,model,v,T,x)
    R = R_GAS
    return log.(x .* (R*T/v)) .+ dAdn
end

function core_logfugacityn(model,v,T,n)
    return core_logfugacityn(model,v,T,n/sum(n))
end

Base.length(model::AbstractHelmholtzModel)=compounds_number(model)

#helmholtz but in total units (volume in m3, n in mol)
function  _helmholtzn(model::AbstractHelmholtzModel,V,T,n)
    sum_n = moles(mol_number(n))
    sum_n< zero(sum_n) && return zero(sum_n)
    x = n*inv(sum_n)
    v = V*inv(sum_n) #this operation seems simple but it propagates the derivative of molar numbers,
    #since core_helmholtz is defined in the sense of molar fractions
    #with that in mind, the chemical potential of a mixture is n*gradient(_helmholtzn) in n
    return sum_n*core_helmholtz(model,v,T,x)
end

function partial_molar_volume(model,v,T,x)
    return partial_molar_property(core_pressure,model,v,T,x)./core_dpdv(model,v,T,x)
end
#helmholtz density, energy in a cubic meter, used in VTC equilibrium
function _helmholtzd(model::AbstractHelmholtzModel,T,C)
    sum_c = sum(C)
    #sum_c < zero(sum_c) && return zero(sum_c)
    return sum_c*core_helmholtz(model,inv(sum_c),T,C/sum_c)
end

#gibbs density, energy in a cubic meter, used in VTC equilibrium
function _gibbsd(model::AbstractHelmholtzModel,T,C)
    sum_c = sum(C)
    sum_c < zero(sum_c) && return zero(sum_c)
    return sum_c*(core_helmholtz(model,inv(sum_c),T,C/sum_c)+core_pressure(model,inv(sum_c),T,C/sum_c))
end

function core_gibbs(model::AbstractHelmholtzModel,v,T,x)
    return core_helmholtz(model,v,T,x)+core_pressure(model,v,T,x)*v
end

function core_gibbs(model::AbstractHelmholtzModel,P,v,T,x)
    return core_helmholtz(model,v,T,x)+P*v
end
#this has the neat property: core_helmholtz(model,v,T,x)= -PV - sum(Ni*μi)
function _chemical_potential(model::AbstractHelmholtzModel,v,T,x)
f = z-> _helmholtzn(model,v,T,z)
return ForwardDiff.gradient(f,x)
end

function _chemical_potentiald(model::AbstractHelmholtzModel,v,T,x)
    f = z-> _helmholtzd(model,T,z/v)
    return ForwardDiff.gradient(f,x)
    end

function core_dpdv(model::AbstractHelmholtzModel,v,T,x)
    return ForwardDiff.derivative(z->core_pressure(model,z,T,x),v)
end

function _d2pdv2(model::AbstractHelmholtzModel,v,T,x)
    return ForwardDiff.derivative(z->_dpdv(model,z,T,x),v)
end

function core_grad_vt(model::AbstractHelmholtzModel,v,T,x)
    d = [v,T]
    f = z-> core_helmholtz(model,z[1],z[2],x)
    return ForwardDiff.gradient(f,d)
end

function core_grad_x(model::AbstractHelmholtzModel,v,T,x)
    f = z-> core_helmholtz(model,v,T,z)
    return ForwardDiff.gradient(f,x)
end

# core_fg functions return value and gradient in one step, reusing calculations
function core_fg_x(model::AbstractHelmholtzModel,v,T,x)
    df = DiffResults.GradientResult(x)
    df = ForwardDiff.gradient!(df,z->core_helmholtz(model,v,T,z),x)
    return (DiffResults.value(df),DiffResults.gradient(df))
end

function core_fg_vt(model::AbstractHelmholtzModel,v,T,x)
    vt = [v,T]
    df = DiffResults.GradientResult(vt)
    df = ForwardDiff.gradient!(df,z->core_helmholtz(model,z[1],z[2],x),vt)
    return (DiffResults.value(df),DiffResults.gradient(df))
end

function core_dfdv(model::AbstractHelmholtzModel,v,T,x)
    return ForwardDiff.derivative(z->core_helmholtz(model,z,T,x),v)
end

function core_dfdt(model::AbstractHelmholtzModel,v,T,x)
    return ForwardDiff.derivative(
        z->core_helmholtz(model,v,z,x),T)
end

function core_hessian_vt(model::AbstractHelmholtzModel,v,T,x)
    f = z->core_helmholtz(model,z[1],z[2],x)
    return ForwardDiff.hessian(f,[v,T])
end

function core_fgh_vt(model::AbstractHelmholtzModel,v,T,x)
    vt = [v T]
    fhelmholtz = z->core_helmholtz(model,z[1],z[2],x)
    df = DiffResults.HessianResult(vt)
    df = ForwardDiff.hessian!(df,fhelmholtz,vt)
    return (DiffResults.value(df),
    DiffResults.gradient(df),
    DiffResults.hessian(df))
end

struct HelmholtzResultData
    value::Array{Float64,1}
    diffresult::DiffResults.DiffResult
end

function _diffdata(model::AbstractHelmholtzModel ,V,T,x,order::Int64 = 2)
    d =  vcat(V,T,x[:])
    #f(d) = core_helmholtz(model,d[1],d[2],d[3:end])
    f = z->  core_helmholtz(model,z[1],z[2],z[3:end])
    res = DiffResults.HessianResult(d)
    if order == 0
    elseif order == 1
    res= ForwardDiff.gradient!(res, f, d)
    elseif order == 2
    res= ForwardDiff.hessian!(res, f, d)
    else
        throw(DomainError(order, """

diffdata can't calculate the derivatives of order $order.
only gradients (order == 1) and hessians (order == 2) are implemented.

"""))
    end
    return HelmholtzResultData(d,res)
end
function diffdata(model::AbstractHelmholtzModel ,v,T,x,order::Int64 = 2)
(v2,T2)=_transformVT(v,T,model.molecular_weight,x)
return _diffdata(model,v2,T2,x,order)
end

_valuevt(df::HelmholtzResultData)  = df.value
_helmholtzvt(df::HelmholtzResultData)  = DiffResults.value(df.diffresult)
core_grad_vt(df::HelmholtzResultData)  = DiffResults.gradient(df.diffresult)[1:2]
core_hessian_vt(df::HelmholtzResultData)  = DiffResults.hessian(df.diffresult)[1:2,1:2]

#
#Pressure, a single derivative is faster
#
function core_pressure(model::AbstractHelmholtzModel,v,T,x)
    return -core_dfdv(model,v,T,x)
return
end

function pressure(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2) = _transformVT(v,T,molecular_weight(model),x)
    return core_pressure(model,v2,T2,x)*1.0u"Pa"
end

function pressure(df::HelmholtzResultData)
    return core_pressure(df)*1.0u"Pa"
end

function core_compressibility_factor(model::AbstractHelmholtzModel,v,T,x)
    return core_pressure(model,v,T,x)*v/(R_GAS*T)
end

function core_compressibility_factor(model::AbstractHelmholtzModel,P,v,T,x)
    return P*v/(R_GAS*T)
end

function compressibility_factor(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2) = _transformVT(v,T,molecular_weight(model),x)
    return core_compressibility_factor(model,v2,T2,x)
end

function compressibility_factor(df::HelmholtzResultData)
    (v,T)=_valuevt(df)
    P= core_pressure(df)
    return P*v/(Unitful.R_GAS*T)
end

#
#Entropy
#
function core_entropy(model::AbstractHelmholtzModel,v,T,x)
    return -core_dfdt(model,v,T,x)
return
end
function core_entropy(df::HelmholtzResultData)
    return -core_grad_vt(df)[2]
end

function entropy(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2) = _transformVT(v,T,molecular_weight(model),x)
    return core_entropy(model,v2,T2,x)*1.0u"J/mol/K"
end

function entropy(df::HelmholtzResultData)
    return core_entropy(df)*1.0u"J/mol/K"
end

#
#enthalpy
#
function core_enthalpy(model::AbstractHelmholtzModel,v,T,x)
    f = core_helmholtz(model,v,T,x)
    df = core_grad_vt(model,v,T,x)
    return f - df[1]*v-df[2]*T
end

function core_enthalpy(df::HelmholtzResultData)
    x = _valuevt(df)
    a = _helmholtzvt(df)
    da= core_grad_vt(df)
    return a - LinearAlgebra.dot(da,x)
end

function enthalpy(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2) = _transformVT(v,T,molecular_weight(model),x)
    return core_enthalpy(model,v2,T2,x)*1.0u"J/mol"
end

function enthalpy(df::HelmholtzResultData)
    return core_enthalpy(df::HelmholtzResultData)*1.0u"J/mol"
end

#
#Internal Energy
#
function core_internal_energy(model::AbstractHelmholtzModel,v,T,x)
    u1 = core_helmholtz(model,v,ForwardDiff.Dual{typeof(T)}(T, one(T)),x)
    return core_helmholtz(model,v,T,x)-core_grad_vt(model,v,T,x)[2]*T
end

function core_internal_energy(df::HelmholtzResultData)
    a = _helmholtzvt(df)
    (v,T)=_valuevt(df)
    da = core_grad_vt(df)
    return a-v*da[2]
end

function internal_energy(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2) = _transformVT(v,T,molecular_weight(model),x)
    return core_internal_energy(model,v2,T2,x)*1.0u"J/mol"
end

function internal_energy(df::HelmholtzResultData)
    return core_internal_energy(df)*1.0u"J/mol"
end

##########
# Cv
##########
function core_isochoric_heat_capacity(model::AbstractHelmholtzModel,v,T,x)
    d2f = core_hessian_vt(model,v,T,x)
    return  -T*d2f[2,2]
    D
end

function core_isochoric_heat_capacity(df::HelmholtzResultData)
    (v,T)= _valuevt(df)
    d2f = core_hessian_vt(df)
    return  -T*d2f[2,2]

end

function isochoric_heat_capacity(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2)=_transformVT(v,T,molecular_weight(model),x)
    return core_isochoric_heat_capacity(model,v2,T2,x)*1.0u"J/mol/K"
end

function isochoric_heat_capacity(df::HelmholtzResultData)
    return core_isochoric_heat_capacity(df)*1.0u"J/mol/K"
end
##########
# Cp,testing
##########
function core_isobaric_heat_capacity(model::AbstractHelmholtzModel,v,T,x)
    d2f = core_hessian_vt(model,v,T,x)
    return  -v*d2f[1,2]-T*d2f[2,2]
end

function core_isobaric_heat_capacity(df::HelmholtzResultData)
    (v,T)=_valuevt(df)
    d2f = core_hessian_vt(df)
    return  -v*d2f[1,2]-T*d2f[2,2]
end

function isobaric_heat_capacity(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2)=_transformVT(v,T,molecular_weight(model),x)
    return     core_isobaric_heat_capacity(model,v2,T2,x)*1.0u"J/mol/K"
end

function isobaric_heat_capacity(df::HelmholtzResultData)
  return core_isobaric_heat_capacity(df::HelmholtzResultData)*1.0u"J/mol/K"
end
#######
# Speed of sound, testing, needs checking
#######
function core_sound_speed(model::AbstractHelmholtzModel,v,T,x)
    d2f = core_hessian_vt(model,v,T,x)
    cp = -v*d2f[1,2]-T*d2f[2,2]
    cv = -T*d2f[2,2]

    m = 0.001*sum(molecular_weight(model).*x)
    return  v*sqrt(d2f[1,1]*(cp/cv)/m)
end

function sound_speed(model::AbstractHelmholtzModel,v,T,x)
    (v2,T2)=_transformVT(v,T,molecular_weight(model),x)
    return     core_sound_speed(model,v2,T2,x)*1.0u"m/s"
end

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

#mixing rules

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

struct GeometricMeanRule <: AbstractMixingRule end
struct ArithmeticMeanRule <: AbstractMixingRule end
struct HarmonicMeanRule <: AbstractMixingRule end
struct PowerMeanRule{T} <: AbstractMixingRule where T<:Real
n::T
end #look at op_exponent below
#
#Abstract mixing rule, generates a mixing rule,based on
# an operation, so mij = xi*xj*op(pi,pj)
#
function mixing_rule(op,x,p)
    N = length(x)
    @boundscheck checkbounds(p,N)
    @inbounds begin
        res1 = zero(eltype(x))
        for i = 1 : N
            res1 += p[i] * x[i]^2
                 for j = 1 : i - 1
                res1 += 2 * x[i] * x[j] * op(p[i], p[j])
            end
        end
    end
    return res1
end
#
#Abstract mixing rule, generates a mixing rule,based on
# an operation, so mij = xi*xj*op(pi,pj)*Aij
#example: mixing_rule(geometric_mean_rule,x,Tc,1.-K)
#
function mixing_rule(op,x,p,A)
    N = length(x)
    checkbounds(A,N,N)
    @boundscheck checkbounds(p,N)
    @inbounds begin
        res1 = zero(eltype(x))
        for i = 1 : N
            res1 += p[i] * x[i]^2
                 for j = 1 : i - 1
                res1 += 2 * x[i] * x[j] * op(p[i], p[j])*A[i,j]
            end
        end
    end
    return res1
end
#
#Abstract asymetric mixing rule
#Adds a Asymetric matrix A_asym, and a op_sim(xi,xj,Aij)
#the mayor example is the GERG2008 equation, where
#op_asym(xi,xj,Aij) = (xi+xj)/(Aij^2*xi + xj)
#
function mixing_rule_asymetric(op,op_asym,x,p,A,A_asym)
    N = length(x)
    checkbounds(A,N,N)
    checkbounds(A_asym,N,N)
    @boundscheck checkbounds(p,N)
    @inbounds begin
        res1 = zero(eltype(x))

        for i = 1 : N
            x[i] !=0 && begin
            res1 += p[i] * x[i]^2
                for j = 1 : i - 1
                    res1 += 2 * x[i] * x[j] * op(p[i], p[j])*A[i,j] *op_asym(x[i],x[j],A_asym[i,j])
                end
            end
        end
        end
    return res1
end

geometric_mean_rule(a,b)=sqrt(a*b)

arithmetic_mean_rule(a,b)=(a+b)/2

harmonic_mean_rule(a,b)=2*a*b/(a+b)

cubic_mean_rule(a,b)=((a^(1/3) + b^(1/3))/2)^3

_power_mean_rule(a,b,n)=((a^(1/n) + b^(1/n))/2)^n

power_mean_rule(n) =(a,b)-> _power_mean_rule(a,b,n)

###############################
#Abstract mixing_matrix
#Creates A symetric matrix from an operation and a vector
#example: A[i,j]=sqrt(T[i]*T[j])
###############################
#
function mixing_matrix(op,p)
    N = length(p)
    A = Array{eltype(p),2}(undef,N,N)
    @boundscheck checkbounds(A,N,N)
    @inbounds begin
        res1 = zero(eltype(p))
        for i = 1 : N
            A[i,i] = p[i]
                for j = 1 : i - 1
                A[i,j] = op(p[i], p[j])
                A[j,i] = op(p[i], p[j])
            end
        end
    end
    return A
end

function mixing_matrix!(op,A,p)
    N = length(p)

    @boundscheck size(A)==(N,N)
    @inbounds begin
        res1 = zero(eltype(p))
        for i = 1 : N
            A[i,i] = p[i]
                for j = 1 : i - 1
                A[i,j] = op(p[i], p[j])
                A[j,i] = op(p[i], p[j])
            end
        end
    end
    return A

end

#mass and molar types for dispatch
#the default is store the fractions and the amount.
#for molar types, the amount is in moles (gmol, as used in the metric system)
#for mass types, the amount is in kg

abstract type AbstractMaterialVector{T} <: AbstractVector{T} end

Base.size(A::AbstractMaterialVector) = size(A.n_values) #size of the vector
Base.length(A::AbstractMaterialVector)=length(A.n_values)
Base.IndexStyle(::Type{AbstractMaterialVector}) = IndexLinear()

struct MolFraction{T} <: AbstractMaterialVector{T}
    n_values::Vector{T}
    n_amount::T
    m_amount::T
    mw::Vector{T}
end

struct MassFraction{T} <: AbstractMaterialVector{T}
    n_values::Vector{T}
    n_amount::T
    m_amount::T
    mw::Vector{T}
end

struct MolNumber{T} <: AbstractMaterialVector{T}
    n_values::Vector{T}
    n_amount::T
    m_amount::T
    mw::Vector{T}
end

struct MassNumber{T} <: AbstractMaterialVector{T}
    n_values::Vector{T}
    n_amount::T
    m_amount::T
    mw::Vector{T}
end

Base.getindex(A::MolNumber, i::Int) =A.n_values[i]*A.n_amount
Base.getindex(A::MolNumber, I::Vararg{Int}) = A.n_values[I...]*A.n_amount
Base.getindex(A::MassNumber, i::Int) = A.n_values[i]*A.mw[i]*A.n_amount*1e-3
Base.getindex(A::MassNumber, I::Vararg{Int}) = A.n_values[I...]*A.mw[I...]*A.n_amount*1e-3
Base.getindex(A::MolFraction, i::Int) =A.n_values[i]
Base.getindex(A::MolFraction, I::Vararg{Int}) =A.n_values[I...]
Base.getindex(A::MassFraction, i::Int) = A.n_values[i]*A.mw[i]*A.n_amount*inv(A.m_amount)*1e-3
Base.getindex(A::MassFraction, I::Vararg{Int}) = A.n_values[I...]*A.mw[I...]*A.n_amount*inv(A.m_amount)*1e-3
core_moles(x::AbstractMaterialVector) = x.n_amount
core_mass(x::AbstractMaterialVector) = x.m_amount
moles(x::AbstractMaterialVector) = x.n_amount*1.0u"mol"
mass(x::AbstractMaterialVector) = x.m_amount*1.0u"kg"
#constructors
function mol_fraction(x::AbstractVector,mw::AbstractVector,n_amount_opt = nothing)
    T = promote_type(eltype(x),eltype(mw))
    n = convert(T,sum(x))
    n_values = convert.(T,x)
    n_values ./=n
    if isnothing(n_amount_opt)
        n_amount = sum(x)
    else
        n_amount = convert(T,n_amount_opt)
    end
    m_amount = LinearAlgebra.dot(x,mw)*1e-3
    return MolFraction{T}(n_values,n_amount,m_amount,convert.(T,mw))
end
function mol_number(x::AbstractVector,mw::AbstractVector,n_amount_opt = nothing)
    T = promote_type(eltype(x),eltype(mw))
    n = convert(T,sum(x))
    n_values = convert.(T,x)
    n_values ./=n
    if isnothing(n_amount_opt)
        n_amount = sum(x)
    else
        n_amount = convert(T,n_amount_opt)
    end
    m_amount = LinearAlgebra.dot(x,mw)*1e-3
    return MolNumber{T}(n_values,n_amount,m_amount,convert.(T,mw))
end

function mass_fraction(x::AbstractVector,mw::AbstractVector,m_amount_opt = nothing)
    T = promote_type(eltype(x),eltype(mw))
    if isnothing(m_amount_opt)
        m_amount = sum(x)
    else
        m_amount = convert(T,m_amount_opt)
    end
    n_values = normalizefrac(x .* inv.(mw))
    n_amount = LinearAlgebra.dot(x,inv.(mw))*1e-3
    return MassFraction{T}(n_values,n_amount,m_amount,convert.(T,mw))
end

function mass_number(x::AbstractVector,mw::AbstractVector,n_amount_opt = nothing)
    T = promote_type(eltype(x),eltype(mw))
    if isnothing(m_amount_opt)
        m_amount = sum(x)
    else
        m_amount = convert(T,m_amount_opt)
    end
    n_values = normalizefrac(x .* inv.(mw) .* 1e3)
    n_amount = LinearAlgebra.dot(x,inv.(mw))*1e3
    return MassNumber{T}(n_values,n_amount,m_amount,convert.(T,mw))
end

mass_number(x::AbstractVector,model::AbstractThermoModel) = mass_number(x,molecular_weight(model))
mass_fraction(x::AbstractVector,model::AbstractThermoModel) = mass_fraction(x,molecular_weight(model))
mol_number(x::AbstractVector,model::AbstractThermoModel) = mol_number(x,molecular_weight(model))
mol_fraction(x::AbstractVector,model::AbstractThermoModel) = mol_fraction(x,molecular_weight(model))

mass_number(model::AbstractThermoModel,x::AbstractVector) = mass_number(x,molecular_weight(model))
mass_fraction(model::AbstractThermoModel,x::AbstractVector) = mass_fraction(x,molecular_weight(model))
mol_number(model::AbstractThermoModel,x::AbstractVector) = mol_number(x,molecular_weight(model))
mol_fraction(model::AbstractThermoModel,x::AbstractVector) = mol_fraction(x,molecular_weight(model))

#conversions
function mass_number(x::AbstractMaterialVector{T}) where T
    return MassNumber{T}(x.n_values,x.n_amount,x.m_amount,x.mw)
end
function mass_fraction(x::AbstractMaterialVector{T}) where T
    return MassFraction{T}(x.n_values,x.n_amount,x.m_amount,x.mw)
end
function mol_fraction(x::AbstractMaterialVector{T}) where T
    return MolFraction{T}(x.n_values,x.n_amount,x.m_amount,x.mw)
end
function mol_number(x::AbstractMaterialVector{T}) where T
    return MolNumber{T}(x.n_values,x.n_amount,x.m_amount,x.mw)
end

####################################
#helmholtz phase
#is used to pass around results in a compact manner,maybe add helmholtz result data here?
#this completely defines a phase in helmholtz energy (v,T)
#temperature in Kelvin, molar volume in m3/mol,and the material is in mol number
struct HelmholtzPhase{T} <: AbstractPhase #this does nothing for now, what functions should implement an abstractPhase?
    volume::T
    temperature::T
    material :: AbstractMaterialVector{T}
end

function helmholtz_phase(model,v,T,x)
TT = promote_type(eltype(x),typeof(v),typeof(T))
v1 = convert(TT,v)
T1 = convert(TT,T)
x1 = mol_number(model,x)
return HelmholtzPhase{TT}(v1,T1,x1)
end

function helmholtz_phase(model,v,T,x::AbstractMaterialVector{_T}) where _T
    TT = promote_type(_T,typeof(v),typeof(T))
    v1 = convert(TT,v)
    T1 = convert(TT,T)
    return HelmholtzPhase{TT}(v1,T1,mol_number(x))
end

core_moles(x) = sum(x) #if a normal vector comes to the solver, this can give you the moles
core_moles(a::HelmholtzPhase)= core_moles(a.material)
core_mass(a::HelmholtzPhase) = core_mass(a.material)
core_mol_volume(a::HelmholtzPhase) = a.volume
core_volume(a::HelmholtzPhase) = a.volume*core_moles(a)
core_mol_density(a::HelmholtzPhase) =inv(a.volume)
core_mass_volume(a::HelmholtzPhase) = a.volume*core_moles(a)/core_mass(a)
core_mass_density(a::HelmholtzPhase) =inv(a.volume*core_moles(a)/core_mass(a))
core_temperature(a::HelmholtzPhase) = a.temperature
core_molecular_weight(a::HelmholtzPhase) = a.material.mw

mol_volume(a::HelmholtzPhase) = core_molar_volume(a)*1.0u"mol/m^3"
mol_density(a::HelmholtzPhase) =core_molar_density(a)*1.0u"m^3/mol"
mass_volume(a::HelmholtzPhase) = core_mass_volume(a)*1.0u"m^3/kg"
mass_density(a::HelmholtzPhase) =core_mass_density(a)*1.0u"kg/m^3"
temperature(a::HelmholtzPhase) = core_temperature(a)*1.0u"K"
molecular_weight(a::HelmholtzPhase) = core_molecular_weight(a).*1.0u"g/mol"
moles(a::HelmholtzPhase)= core_moles(a)*1.0u"mol"
mass(a::HelmholtzPhase) = core_mass(a)*1.0u"kg"

mol_fraction(a::HelmholtzPhase) = mol_fraction(a.material)
mol_number(a::HelmholtzPhase) = mol_number(a.material)
mass_fraction(a::HelmholtzPhase) = mass_fraction(a.material)
mass_number(a::HelmholtzPhase) = mass_number(a.material)
_unpack_helmholtz(a::HelmholtzPhase) = (core_molar_volume(a),core_temperature(a),mol_fraction(a))

for op = (:_helmholtzn,    :core_helmholtz, :core_gibbs,    :core_grad_vt,    :core_grad_x,
    :core_fg_x,    :core_fg_vt,    :core_dfdv,    :core_dfdt,    :core_hessian_vt,
    :core_fgh_vt,    :core_pressure,    :compressibility_factor,    :core_entropy,
    :core_enthalpy,    :core_internal_energy,    :core_isochoric_heat_capacity,
    :core_isobaric_heat_capacity,    :core_sound_speed,
    :pressure,     :entropy,    :enthalpy,    :internal_energy,    :isochoric_heat_capacity,
    :isobaric_heat_capacity,    :sound_speed, :core_dpdv)
    @eval $op(model::AbstractHelmholtzModel,a::HelmholtzPhase) = $op(model,core_mol_volume(a),core_temperature(a),mol_fraction(a))
    @eval $op(model::AbstractHelmholtzModel,v,T,x::Union{MassFraction,MassFraction,MassNumber}) = $op(model,v,T,mol_fraction(x))
end


text:
function position_mod(x::Float64, L::Float64)
	return mod(x, L) # For characterization we need to do full modulus operation.
end

ext: jl


text:
using ALASMPCCmod
using MPCCmod
using Relaxation
using MPCC2DPlot

r=1.0;s=1.0;t=1.0;
ALASMPCCmod_success=true
println("On teste le package ALASMPCC")

nlp=ADNLPModel(f,[1.0;-1.0], lvar=-Inf*ones(2), uvar=Inf*ones(2), c=c,y0=ones(1),lcon=lcon,ucon=Inf*ones(1))
mpcc = MPCCmod.MPCC(nlp,G,H,1)
alas = ALASMPCCmod.ALASMPCC(mpcc,r,s,t)

#SlackComplementarityProjection(x::Vector,nb_comp::Int64,r::Float64,s::Float64,t::Float64)
#y1=ALASMPCC.SlackComplementarityProjection([1.0;-1.0;-2.0;0.5],1,r,s,t)
#y2=ALASMPCC.SlackComplementarityProjection([1.0;-1.0;1.0;-3.0],1,r,s,t)
#y3=ALASMPCC.SlackComplementarityProjection([1.0;-1.0;2.0;3.0],1,r,s,t)
#if !(y1[3]>=-r && y2[4]>=-r && y3[4]==Relaxationmod.psi(y3[3],r,s,t))
# println(y1,y2)
#
# ALASMPCCmod_success=false
#end

r=1.0;s=1.0;t=0.0;
nlp=ADNLPModel(f,[2.0;-1.0], lvar=-Inf*ones(2), uvar=Inf*ones(2), c=c,y0=ones(1),lcon=lcon,ucon=Inf*ones(1))
mpcc = MPCCmod.MPCC(nlp,G,H,1)
alas = ALASMPCCmod.ALASMPCC(mpcc,r,s,t)

println("Exemple 1")
#output=ALASMPCC.solveALAS(mpcc,r,s,t)
output=ALASMPCCmod.solveALAS(alas)
println("solution :", output[1])
println("multiplicateurs :", output[2]," ",output[3]," ",output[4])
println("nb itérations :", length(output[5])/4)
nb_ite=Int(length(output[5])/4)
x1=zeros(nb_ite);x2=zeros(nb_ite)
for i=1:nb_ite
 x1[i]=output[5][4*(i-1)+1] ; x2[i]=output[5][4*(i-1)+2]
end
 PyPlot.plot(x1, x2, color="red", linewidth=2.0, linestyle="-")
 PyPlot.xlabel("x1")
 PyPlot.ylabel("x2")
 PyPlot.title("min x1-x2 s.t. x2<=1,x1>=-1,x2>=-1,(x1-1)=0 ou (x2-1)=0")

#println("Exemple 2:")
r=1.0;s=1.0;t=0.0;
nlp=ADNLPModel(f,[1.0;1.5], lvar=-Inf*ones(2), uvar=Inf*ones(2), c=c,y0=ones(1),lcon=lcon,ucon=Inf*ones(1))
mpcc = MPCCmod.MPCC(nlp,G,H,1)
alas = ALASMPCCmod.ALASMPCC(mpcc,r,s,t)
#output=ALASMPCC.solveALAS(mpcc,r,s,t)

# Un exemple simple pour tester :
#println("Exemple 3:")
# f2(x)=x[1]^2-x[1]*x[2]+1.0/3.0*x[2]^2-2*x[1]
# G2(x)=x[1]
# H2(x)=x[2]
# cc(x)=[x[1]+10]
# lcon2=zeros(1)
r=1.0;s=1.0;t=0.0;
#nlp2=ADNLPModel(f2,[-1.0;1.0], lvar=-Inf*ones(2), uvar=Inf*ones(2), c=cc,y0=ones(1),lcon=lcon2,ucon=Inf*ones(1))
#mpcc2 = MPCCmod.MPCC(nlp2,G2,H2,1)
#alas = ALASMPCC.ALASMPCC(mpcc,r,s,t)
#output=ALASMPCC.solveALAS(mpcc2,r,s,t)

#MPCC2DPlot.RelaxationPlot([1.0 -1.0;1.0 1.5],r,s,t,-2.0,2.0)

#Bilan : donner une sortie bilan
if ALASMPCCmod_success==true
 println("ALASMPCC.jl passes the test !")
else
 println("ALASMPCC.jl contains some error")
end


text:
using BinaryBuilder

# Collection of sources required to build WABT
sources = [
    "https://github.com/WebAssembly/wabt/archive/1.0.0.tar.gz" =>
    "a5d4cfb410fbe94814ed8ae67a2c356c4ea39d26578ca5b48a8d7ede2a0e08eb",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd wabt-1.0.0/
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_BUILD_TYPE=gcc-release -DBUILD_TESTS=OFF -DCMAKE_CXX_FLAGS="-static -no-pie" -DCMAKE_C_FLAGS="-static -no-pie"
else
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_BUILD_TYPE=gcc-release -DBUILD_TESTS=OFF
fi

make -j8
make install

if [ $target = "x86_64-w64-mingw32" ]; then
  echo hello
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Windows(:x86_64),
    BinaryProvider.Windows(:i686),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:i686, :glibc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.Linux(:powerpc64le, :glibc, :blank_abi),
    BinaryProvider.Linux(:i686, :musl, :blank_abi),
    BinaryProvider.Linux(:x86_64, :musl, :blank_abi),
    BinaryProvider.Linux(:aarch64, :musl, :blank_abi),
    BinaryProvider.Linux(:armv7l, :musl, :eabihf)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "", :wast_sugar),
    ExecutableProduct(prefix, "", :wasm_objdump),
    ExecutableProduct(prefix, "", :wasm_interp),
    ExecutableProduct(prefix, "", :wasm_link),
    ExecutableProduct(prefix, "", :wasm_opcodecnt),
    ExecutableProduct(prefix, "", :wast2wasm),
    ExecutableProduct(prefix, "", :wasm2wast)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

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

max_stars_repo_name: MikeInnes/WABTBuilder

max_issues_repo_name: MikeInnes/WABTBuilder

max_forks_repo_name: MikeInnes/WABTBuilder


text:
# p16.jl - Poisson eq. on [-1,1]x[-1,1] with u=0 on boundary

# Set up grids and tensor product Laplacian and solve for u:
N = 24; (D,x) = cheb(N); y = x;
xx = x[2:N]; yy = y[2:N];
f = 10*sin.(8*xx'.*(yy-1));
D2 = D^2; D2 = D2[2:N,2:N]; I = eye(N-1);
L = kron(I,D2) + kron(D2,I);                       # Laplacian
figure(1); clf(); spy(L);
tic(); u = L\f[:]; toc();          # solve problem and watch the clock

# Reshape long 1D results onto 2D grid (flipping orientation):
uu = zeros(N+1,N+1); uu[N:-1:2,N:-1:2] = reshape(u,N-1,N-1);
value = uu[Int(3N/4+1),Int(3N/4+1)];

# Interpolate to finer grid and plot:
xxx = yyy = -1:.04:1;
s = interpolate((x[end:-1:1],y[end:-1:1]),uu,Gridded(Linear()));
uuu = s[xxx,yyy];
figure(2); clf(); surf(xxx,yyy,uuu,rstride=1,cstride=1);
xlabel("x"); ylabel("y"); zlabel("u"); view(-37.5,30);
text3D(.4,-.3,-.3,"\$u(2^{-1/2},2^{-1/2})\$ = $(signif(value,11))",fontsize=9);


text:
using Pkg
Pkg.activate(pwd())

# # Arrays
# ## Vectors

v = [1, 2, 3, 4, 5, 6, 7, 8] # or equivalently v = [1; 2; 3; 4; ...]

#+

typeof(v)
ndims(v)
eltype(v)
size(v)
length(v)

#+

v[3]
v[begin] # the first element
v[end] # the last element
v[[2, 3]]

#+

range(1; stop = 10, step = 2) # or equivalently range(1, 10; step = 2)
1:2:10

#+

v[1:3] # the first three elements
v[1:2:end] # select all elements with odd index
v[:] # all elements

#+

v = [1,2,3]

append!(v, 4)
append!(v, [5,6])
append!(v, 7:8)
append!(v, 3.0)
append!(v, 3.1415)

#+

isinteger(3.0)

#+

v = Float64[1, 2, 3]
append!(v, 3.1415)

#+

v = [1, 2, 3, 4]
v[2] = 4
v

#+

v[3:4] .= 11
v

# ### Exercise:
# Create a vector of positive integers that contains all odd numbers smaller than `10`.
# Then change the first element to `4` and the last two elements to `1`.
#
# ---
# ### Solution:

v = [1,3,5,7,9]

collect(1:2:9)
Vector(1:2:9)

v[1] = 4
v[end-1:end] .= 1
v

# ---
#
# ## Matrices

m = [1  2  3  4; 5  6  7  8]

typeof(m)
eltype(m)
ndims(m)
size(m)
length(m)

#+

m[1] # the first element, equivalent to m[begin]
m[2] # the second element
m[end-1] # the last element
m[1, 2]
m[1, [2, 3]] # the second and third element in the first row
m[1:3] # the first three elements according to linear indexing
m[:, 1:3] # the first three columns
m[1, :] # the first row
m[:] # all elements

#+

hcat(m, m)
vcat(m, m)

#+

cat(m, m; dims = 2) # equivalent to hcat(m, m)
cat(m, m; dims = 1) # equivalent to vcat(m, m)

#+

v = [11, 12]
vcat(m, v)

# ### Exercise:
# Create two vectors: vector of all odd positive integers smaller than `10` and vector of
# all even positive integers smaller than `10`. Then concatenate these two vectors
# horizontally and fill the third row with `4`.
#
# ---
# ### Solution:

v1 = collect(1:2:9)
v2 = collect(2:2:10)

M = hcat(v1, v2)
M[3,:] .= 4
M

# ---
#
# ## `N`-dimensional arrays

A = zeros(3, 5, 2) # equivalent to A = zeros((3, 5, 2))
B = zeros(Int64, 3, 5, 2)  # equivalent to B = zeros(Int64, (3, 5, 2))

typeof(A)
eltype(A)
ndims(A)
size(A)
length(A)

#+

B[1] = 1 # assign 1 to the first element
B[1, 2, 2] = 2 # assign 2 to the element at position (1,2,2)
B[2,:,1] .= 4

B

#+

ones(Float32, 2, 3, 1)
fill(1.234, 2, 2)

# ### Exercise:
# Create three matrices with the following properties:
# - Matrix `A` is of size `2x3`, and all its elements equal 0.
# - Matrix `B` is of size `2x3x1`, and all its elements equal 1.
# - Matrix `C` is of size `2x3`, and all its elements equal 2.
# Concatenate these three matrices along the third dimension.
#
# **Hint:** use the `cat` function and the keyword `dims`.
#
# ---
# ### Solution:

A = zeros(2, 3)
B = ones(2, 3, 1)
C = fill(2, 2, 3)

cat(A, B, C; dims = 3)

# ---
#
# ## Broadcasting

a = [1,2,3] # column vector

a .-= 4 # from each element of vector subtracts 4
a -= 1

#+

abs.(a)
sum(exp.(sqrt.(abs.(a .- 1)))./2)

#+

a = [1,2,3] # column vector
b = [4,5,6] # column vector

a * b
a' * b
a * b'
a .* b

# ### Exercise:
# Construct a matrix whose elements are given by the following formula
# $$
# A_{i, j} = \frac{1}{2}\exp\{(B_{i, j} + 1)^2\}, \quad i \in \{1, 2\}, \; j \in
# \{1, 2, 3\}
# $$
# where the matrix $B$ is defined by

B = [
    -1  0  2;
    2  -3  1;
]

#
# ---
# ### Solution:

A = exp.((B .+ 1) .^ 2) ./ 2

A = @. exp((B + 1) ^ 2) / 2

A = zeros(2, 3);
for i in 1:length(A)
    A[i] = exp((B[i] + 1)^2)/2
end
A

# ---
#
# ## Views

A = [1 2 3; 4 5 6]
B = A

B[2] = 42

A == B

#+

C = copy(A)
C[4] = 10

A == C

#+

D = A[1:2, 1:2]
D[1] = 15

D == A[1:2, 1:2]

#+

E = view(A, 1:2, 1:2)
E = @view A[1:2, 1:2]

E[4] = 78

E == A[1:2, 1:2]

#+

typeof(E)

#+

A = [1 2 3; 4 5 6]
A_view = @view A[:, :]

sum(A)
sum(A_view)
minimum(A; dims = 1)
minimum(A_view; dims = 1)

#+

A = [1 2 3; 4 5 6];

sum(exp.(sqrt.(abs.(@view(A[1, :]) .- @view(A[2, :]))))./2)
@views sum(exp.(sqrt.(abs.(A[1, :] .- A[2, :])))./2)


text:
##--------------------------------##
##编写求解器函数
##1、SVD-l1
##--------------------------------##
#找到占总比重大于weight前几个特征值,svd分解后s本身已经是从大到小排序完的
export solver

function find_num_max(s, weights)
    @assert typeof(s) <: Vector
    ss = abs.(s)
    weight = 0.0
    total = sum(ss)
    j = 0::Int
    for i in ss
       weight += i/total
       j += 1
       weight >= weights && break
    end
    return j
end

##求解Postion结构体的DOA
function solver(x::Position, fast)
    ex = Expr(:call, x.solver, x, fast)
    eval(ex)
end
##调用前先将data关于x的元素全部清0
function svd_l1!(data::Position,fast = false; weight = 0.95, hyperparameter = 0.8, )

    @assert size(data.real_A, 2) == size(data.real_x, 1)
    @assert size(data.real_A, 2) == size(data.imag_x, 1)

    real_A = view(data.real_A, :, :)
    imag_A = view(data.iamg_A, :, :)
    y = data.y

    f = svd(y)
    num_K = find_num_max(f.S, weight)
    IK = zeros(eltype(real_A), num_Ksize(real_A, 1), num_K)
    Threads.@threads for i in 1:num_K
        IK[i, i] = one(eltype(IK))
    end
    Ysv = y* f.V* IK
    real_Ysv = real.(Ysv)
    imag_Ysv = imag.(Ysv)
    n = size(real_A, 2)

    if fast
       ##采用局部细化网格加速计算，但我也不知道有没有问题hh
       ##在原本划分网格的程度上每隔十个点下采样，然后再是原始网格
       seq = Array{Bool, 1}(undef, n)
       fill!(seq, false)
        for i in 1:2
            if i == 1
                real_A1 = view(real_A, :, 1:10:n)
                imag_A1 = view(imag_A, :, 1:10:n)
            elseif i == 2
                real_A1 = view(real_A, :, seq)
                imag_A1 = view(imag_A, :, seq)
            end
            n1 = size(real_A1, 2)

            model = JuMP.Model(COSMO.Optimizer)
            @variable(model, s_real[1:n1, 1:num_K])
            @variable(model, s_imag[1:n1, 1:num_K])
            @variable(model, 0 <= p)
            @variable(model, z_real[1:length(Ysv)])
            @variable(model, z_imag[1:length(Ysv)])
            @constraint(model, z_real .== (real_Ysv - real_A1*s_real + imag_A1*s_imag)[1:end])
            @constraint(model, z_imag .== (imag_Ysv - imag_A1*s_real - real_A1*s_imag)[1:end])
            @constraint(model, [p; z_real; z_imag] in SecondOrderCone())

            @variable(model, r[1:n1])
            for i in 1:n1 #
                @constraint(model, [r[i]; s_real[i, :]; s_imag[i, :]] in SecondOrderCone())
            end
            @variable(model, 0 <= q)
            @constraint(model, q >= sum(r))
            @objective(model, Min, p^2+hyperparameter*q)
            JuMP.optimize!(model)

            if i == 1
                let
                    Ssv = sum(sqrt.((value.(s_real)).^2 + (value.(s_imag)).^2), dims = 2)
                    power = 10*log10.(Ssv ./ findmax(Ssv)[1]) ##功率大于-10的前几个点的周围
                    index = partialsortperm(power[:, 1], 1:(num_K+1), rev = true)
                    @views Threads.@threads for i in 1:length(index)
                        max = index[i]+5 > n ? n : index[i] + 5
                        min = index[i]-5 < 1 ? 1 : index[i] - 5
                        @. seq[min:1:max] = true
                    end
                end
            elseif i == 2
                ##data中x的大小肯定大于细化网格结果大小
                ##将计算结果置入data中，重新计算seq
                @. data.real_x[seq, 1:num_K] = value(s_real)
                @. data.imag_x[seq, 1:num_K] = value(s_imag)
                let
                    Ssv = sum(sqrt.((value.(s_real)).^2 + (value.(s_imag)).^2), dims = 2)
                    power = 10*log10.(Ssv ./ findmax(Ssv)[1]) ##功率大于-10的前几个点的周围
                    index = partialsortperm(power[:, 1], 1:num_K, rev = true)
                    fill!(seq, false)
                    @. seq[index] = true
                end

            else
                error("i should be 1 or 2")
            end

        end
        return seq #返回需要的序列用于绘图
    else
         #那这每次运行模型都是重新，吃完晚饭回来去JuMP文档中查看更新参数的功能

         model = JuMP.Model(COSMO.Optimizer)
         @variable(model, s_real[1:n, 1:num_K])
         @variable(model, s_imag[1:n, 1:num_K])
         @variable(model, 0 <= p)
         @variable(model, z_real[1:length(Ysv)])
         @variable(model, z_imag[1:length(Ysv)])
         @constraint(model, z_real .== (real_Ysv - real_A*s_real + imag_A*s_imag)[1:end])
         @constraint(model, z_imag .== (imag_Ysv - imag_A*s_real - real_A*s_imag)[1:end])

         @constraint(model, [p; z_real; z_imag] in SecondOrderCone())

         @variable(model, r[1:n])
         for i in 1:n
            @constraint(model, [r[i]; s_real[i, :]; s_imag[i, :]] in SecondOrderCone())
         end

         @variable(model, 0<=q)
         @constraint(model, q >= sum(r))
         @objective(model, Min, p^2+hyperparameter*q)

         JuMP.optimize!(model)

        data.real_x[:, 1:num_K] = value.(s_real)
        data.imag_x[:, 1:num_K] = value.(s_imag) #将结果写入
        #  value.(s_imag)
    end

    return true
end


text:
@testset "entropy.jl" begin

    Random.seed!(6)

    @testset "entropy" begin

        @test entropy([0.5, 0.5]) ≈ 1
        @test entropy([0.25, 0.25, 0.25, 0.25]) ≈ 2

        pvec = rand(10)
        pvec ./= sum(pvec)
        ρ = Diagonal(pvec)

        @test entropy(ρ) ≈ entropy(Matrix(ρ))
        @test entropy(Hermitian(Matrix(ρ))) ≈ entropy!(Matrix(ρ))
    end

    @testset "entanglement_entropy" begin
        bell_state = normalize(up⊗down + down⊗up)
        @test entanglement_entropy(up⊗up) ≈ 0
        @test entanglement_entropy(bell_state) ≈ 1
        @test entanglement_entropy(bell_state ⊗ down) ≈ 0
        @test entanglement_entropy(bell_state ⊗ down, 2) ≈ 1
        @test entanglement_entropy(bell_state ⊗ bell_state, 1) ≈ 1
        @test isapprox(entanglement_entropy(bell_state ⊗ bell_state, 2), 0; atol=1e-15)
        @test entanglement_entropy(bell_state ⊗ bell_state, 3) ≈ 1
    end
end


text:
mutable struct VarIsing
    uin::Vector{Float64}
    uout::Vector{Ptr{Float64}}
    tJ::Vector{Float64}
    H::Float64
    htot::Float64
end

VarIsing() = VarIsing(Vector{Float64}(), Vector{Ptr{Float64}}(), Vector{Float64}(), 0, 0)

mag(v::VarIsing) = tanh(v.htot)
deg(v::VarIsing) = length(v.uin)

function Base.show(io::IO, v::VarIsing)
    print(io, "VarIsing(deg=$(deg(v)), H=$(v.H))")
end

mutable struct FGIsing
    N::Int
    vnodes::Vector{VarIsing}
    adjlist::Vector{Vector{Int}}
    mags::Vector{Float64}
end

function FGIsing(net::Network; T=1)
    @assert has_eprop(net, "J") || has_gprop(net, "J")
    @assert has_vprop(net, "H") || has_gprop(net, "H")

    hasconstJ = !has_eprop(net, "J")
    hasconstH = !has_vprop(net, "H")

    adjlist = adjacency_list(net)
    N = nv(net)
    vnodes = [VarIsing() for i=1:N]
    if hasconstJ
        J = gprop(net, "J") / T
    else
        J = [[eprop(net, e)["J"] / T for e in edges(net, i)] for i=1:N]
    end
    for (i, v) in enumerate(vnodes)
        resize!(v.uin, length(adjlist[i]))
        resize!(v.uout, length(adjlist[i]))
        resize!(v.tJ, length(adjlist[i]))
    end

    for (i, v) in enumerate(vnodes)
        for (ki, j) in enumerate(adjlist[i])
            v.uin[ki] = 0
            kj = findfirst(==(i), adjlist[j])
            vnodes[j].uout[kj] = getref(v.uin, ki)
            if hasconstJ
                v.tJ[ki] = tanh(J)
            else
                v.tJ[ki] = tanh(J[i][ki])
            end
        end
        if hasconstH
            v.H = gprop(net, "H") / T
        else
            v.H = vprop(net, i)["H"] / T
        end
        v.htot = sum(v.uin) + v.H
    end

    mags = mag.(vnodes)

    FGIsing(N, vnodes, adjlist, mags)
end

function initrand!(g::FGIsing; μ=0, σ=1)
    for v in g.vnodes
        for k=1:deg(v)
            v.uin[k] = μ + σ * randn()
        end
        v.htot = sum(v.uin) + v.H
    end
end

function getmess(g::FGIsing, i, j)
    ki = findfirst(==(j), g.adjlist[i])
    return g.vnodes[i].uout[ki][]
end

function update!(v::VarIsing)
    @unpack uin, uout, tJ = v
    v.htot = sum(uin) + v.H
    for k=1:deg(v)
        hcav = v.htot - uin[k]
        ucav = atanh(tJ[k] * tanh(hcav))
        uout[k][] = ucav
    end
end

function oneBPiter!(g::FGIsing)
    for i in randperm(g.N)
        update!(g.vnodes[i])
    end
    Δ = 0.
    for i=1:g.N
        m = mag(g.vnodes[i])
        Δ = max(Δ, abs(g.mags[i] - m))
        g.mags[i] = m
    end
    return Δ
end

function converge!(g::FGIsing; maxiters=1000, ϵ=1e-6, verbose=true)
    it = 0
    while it < maxiters
        it += 1
        Δ = oneBPiter!(g)
        verbose && @printf("it=%d  Δ=%.2g \n", it, Δ)
        if Δ < ϵ
            verbose && println("Converged!")
            break
        end
    end
    return it
end

function corr_conn_nn(g::FGIsing, i::Int, j::Int)
    @unpack N, vnodes, adjlist = g
    vi = vnodes[i]
    vj = vnodes[j]
    ki = findfirst(==(j), adjlist[i])
    kj = findfirst(==(i), adjlist[j])
    @assert(ki >0)
    @assert(kj >0)

    tJ = vi.tJ[ki]
    mij = tanh(htot(vi) - vi.uin[ki])
    mji = tanh(htot(vj) - vj.uin[kj])

    c = tJ * (1-mij^2)*(1-mji^2) / (1+tJ * mij * mji)
    return c
end

corr_disc_nn(g::FGIsing,i::Int,j::Int) = corr_conn_nn(g,i,j) + mag(g.vnodes[i])*mag(g.vnodes[j])

function corr_disc_nn(g::FGIsing)
    corrs = [zeros(length(g.adjlist[i])) for i=1:g.N]

    for i=1:g.N
        for (ki, j) in enumerate(g.adjlist[i])
            corrs[i][ki] = corr_disc_nn(g, i, j)
        end
    end
    return corrs
end

function run_bp(net::Network;
                maxiters=1000, #max bp iters
                ϵ=1e-6, # stopping criterium
                T=1,    # temperature
                μ=0,    # mean init messages
                σ=1,    # std init messages
                verbose=true,
                )
    g = FGIsing(net; T)
    initrand!(g; μ, σ)
    iters = converge!(g; maxiters, ϵ, verbose)

    res = (bpgraph = g,
           mags = g.mags,
           iters = iters)
    return res
end


text:
#!/home/vagrant/.julia/bin/julia

println("Julia: installing IJulia kernel...")

Pkg.add("IJulia")

max_stars_repo_name: infiRD/Scibox

max_issues_repo_name: infiRD/Scibox

max_forks_repo_name: infiRD/Scibox


text:
# SPDX-License-Identifier: MIT

const D = 5 # vector/matrix dimension
const B = 2 # batch size
const S = 3 # number of states
const K = 3 # Number of emissions
const N = 7 # number of frames

const T = Float32
const SF = LogSemiring{T}

function makefsm(SF, S)
    fsm = VectorFSM{SF}()

    prev = addstate!(fsm, 1, initweight = one(SF))
    addarc!(fsm, prev, prev)
    for s in 2:S
        fw = s == S ? one(SF) : zero(SF)
        state = addstate!(fsm, s, finalweight = fw)
        addarc!(fsm, prev, state)
        addarc!(fsm, state, state)
        prev = state
    end
    fsm |> renormalize
end

function forward(A, init, lhs::AbstractMatrix{T}) where T
    N = size(lhs, 2)
    log_α = similar(lhs)
    log_α[:, 1] = lhs[:, 1] + init
    for n in 2:N
        log_α[:, n] = lhs[:,n] + dropdims(logsumexp(A .+ log_α[:,n-1], dims = 1),
                                          dims = 1)
    end
    log_α
end

function backward(Aᵀ, final, lhs::AbstractMatrix{T}) where T
    N = size(lhs, 2)
    log_β = similar(lhs)
    log_β[:, end] = final
    for n in N-1:-1:1
        log_β[:, n] = dropdims(logsumexp(Aᵀ .+ log_β[:,n+1] .+ lhs[:,n+1], dims = 1),
                               dims = 1)
    end
    log_β
end

function forward_backward(A, Aᵀ, init, lhs)
    final = A[1:end-1,end]

    # ignore the final state
    A = A[1:end-1,1:end-1]
    Aᵀ = Aᵀ[1:end-1,1:end-1]
    init = init[1:end-1]

    log_α = forward(A, init, lhs)
    log_β = backward(Aᵀ, final, lhs)
    log_γ = log_α .+ log_β
    sums = logsumexp(log_γ, dims = 1)
    exp.(log_γ .- sums), minimum(sums)
end

if CUDA.functional()
    @testset "linalg" begin
        D = 5

        x = ones(T, D)
        y = spzeros(T, D)
        y[2] = T(3)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        out = similar(xg)
        MarkovModels.Inference.elmul_svdv!(out, yg, xg)
        @test all(val.(out) .≈ x .* y)

        x = ones(T, D)
        y = spzeros(T, D)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        out = similar(xg)
        MarkovModels.Inference.elmul_svdv!(out, yg, xg)
        @test all(val.(out) .≈ x .* y)

        for Yᵀ in [spdiagm(y), spzeros(T, length(y))]
            Yᵀ = spdiagm(y)
            Yᵀg = CuSparseMatrixCSC(Yᵀ)
            Yg = CuSparseMatrixCSR( Yᵀg.colPtr, Yᵀg.rowVal, Yᵀg.nzVal, (D,D))
            out = similar(xg)
            MarkovModels.Inference.mul_smdv!(out, Yg, xg)
            @test all(val.(out) .≈ transpose(Yᵀ) * x)
        end

        x = ones(T, D)
        y = spzeros(T, D)
        X = ones(T, D, D)
        Y = spzeros(T, D, D)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        Xg = CuArray(X)
        Yg = CuSparseMatrixCSC(Y)
        Yg = CuSparseMatrixCSR(Yg.colPtr, Yg.rowVal, Yg.nzVal, Yg.dims)
        out = similar(Xg)
        MarkovModels.Inference.elmul_svdm!(out, yg, Xg)
        @test all(val.(out) .≈ X .* reshape(y, :, 1))

        x = ones(T, D)
        y = spzeros(T, D)
        X = ones(T, D, D)
        Y = spzeros(T, D, D)
        Y[1,1] = T(1)
        Y[1,2] = T(2)
        Y[2,1] = T(3)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        Xg = CuArray(X)
        Yg = CuSparseMatrixCSC(Y)
        Yg = CuSparseMatrixCSR(Yg.colPtr, Yg.rowVal, Yg.nzVal, Yg.dims)
        out = similar(Xg)
        MarkovModels.Inference.elmul_svdm!(out, yg, Xg)
        @test all(val.(out) .≈ X .* reshape(y, :, 1))

        A = spzeros(3, 3)
        Ag = CuSparseMatrixCSC(A)
        Ag = CuSparseMatrixCSR(Ag.colPtr, Ag.rowVal, Ag.nzVal, Ag.dims)
        Y = ones(3, 3)
        Yg = CuArray(Y)
        out = similar(Yg)
        MarkovModels.Inference.mul_smdm!(out, Ag, Yg)
        @test all(val.(out) .≈ A' * Y)

        A = spzeros(3, 3)
        A[1, 1] = 1/2
        A[1, 2] = 1/2
        A[2, 2] = 1/2
        A[2, 3] = 1/2
        A[3, 3] = 1/2
        Ag = CuSparseMatrixCSC(A)
        Ag = CuSparseMatrixCSR(Ag.colPtr, Ag.rowVal, Ag.nzVal, Ag.dims)
        Y = ones(3, 3)
        Yg = CuArray(Y)
        out = similar(Yg)
        MarkovModels.Inference.mul_smdm!(out, Ag, Yg)
        @test all(val.(out) .≈ A' * Y)

        A = ones(T, 3, 2)
        Ag = CuArray(A)
        B = ones(T, 2, 3)
        Bg = CuArray(B)
        out = similar(Ag, 3, 3)
        MarkovModels.Inference.mul_dmdm!(out, Ag, Bg)
        @test all(val.(out) .≈ A * B)
    end

    @testset "array" begin
        X₁ = spzeros(T, 4, 3)
        X₂ = spzeros(T, 3, 2)

        X₁[1, 1] = T(2)
        X₁[3, 3] = T(3)
        X₁[4, 2] = T(-1)

        X₂[2, 1] = T(2)
        X₂[2, 2] = T(3)

        Y = blockdiag(X₁, X₂)

        gX₁ = CuSparseMatrixCSR(cu(X₁))
        gX₂ = CuSparseMatrixCSR(cu(X₂))
        gY = blockdiag(gX₁, gX₂)

        @test all(val.(gY) .≈ val.(Y))

        x₁ = spzeros(T, 3)
        x₂ = spzeros(T, 2)

        x₁[2] = T(1)
        x₂[1] = T(2)
        y = vcat(x₁, x₂)

        gx₁ = cu(x₁)
        gx₂ = cu(x₂)
        gy = vcat(gx₁, gx₂)

        gy = CUDA.@allowscalar Array(gy)
        @test all(gy .≈ y)
    end
end

@testset "forward_backward" begin
    lhs = ones(T, S, N)
    label2pdfid = Dict(1 => 1, 2 => 2, 3 => 3)

    fsm = makefsm(SF, S)
    mfsm = MatrixFSM(fsm, label2pdfid)

    γ_ref, ttl_ref = forward_backward(
        val.(mfsm.T),
        val.(mfsm.Tᵀ),
        val.(mfsm.π),
        lhs
    )

    γ_scpu, ttl_scpu = pdfposteriors(mfsm, lhs)
    @test all(γ_ref .≈ γ_scpu)
    @test ttl_ref ≈ ttl_scpu

    if CUDA.functional()
        mfsm = mfsm |> gpu
        γ_sgpu, ttl_sgpu = pdfposteriors(mfsm, CuArray(lhs))
        @test all(γ_ref .≈ val.(γ_sgpu))
        @test ttl_ref ≈ ttl_sgpu
    end
end

@testset "batch forward_backward" begin
    lhs = ones(T, S, 7, 2)
    seqlengths = [5, 7]

    label2pdfid = Dict(1 => 1, 2 => 2, 3 => 3)

    fsm = makefsm(SF, S)
    mfsm = MatrixFSM(fsm, label2pdfid)

    γ_ref1, ttl_ref1 = forward_backward(
        val.(mfsm.T),
        val.(mfsm.Tᵀ),
        val.(mfsm.π),
        lhs[:, 1:5, 1]
    )

    γ_ref2, ttl_ref2 = forward_backward(
        val.(mfsm.T),
        val.(mfsm.Tᵀ),
        val.(mfsm.π),
        lhs[:, 1:7, 1]
    )

    mfsms = union([mfsm for i in 1:B]...)
    γ_scpu, ttl_scpu = pdfposteriors(mfsms, lhs, seqlengths)
    @test all(γ_ref1 .≈ γ_scpu[:,1:5,1])
    @test ttl_ref1 ≈ ttl_scpu[1]
    @test all(γ_ref2 .≈ γ_scpu[:,1:7,2])
    @test ttl_ref2 ≈ ttl_scpu[2]
    @test all(γ_scpu[:,6:7,1] .== zero(T))

    if CUDA.functional()
        mfsm = mfsm |> gpu
        mfsms = union([mfsm for i in 1:B]...)
        γ_sgpu, ttl_sgpu = pdfposteriors(mfsms, CuArray(lhs), seqlengths)
        @test all(γ_ref1 .≈ val.(γ_sgpu)[:,1:5,1])
        @test ttl_ref1 ≈ val.(ttl_sgpu)[1]
        @test all(γ_ref2 .≈ val.(γ_sgpu)[:,1:7,2])
        @test ttl_ref2 ≈ val.(ttl_sgpu)[2]
        @test all(val.(γ_sgpu)[:,6:7,1] .== zero(T))
    end
end

@testset "bestpath" begin
    fsm = VectorFSM{SF}()

    s1 = addstate!(fsm, "a", initweight = one(SF))
    s2 = addstate!(fsm, "b")
    s3 = addstate!(fsm, "c")
    s4 = addstate!(fsm, "d", finalweight = one(SF))

    addarc!(fsm, s1, s2)
    addarc!(fsm, s2, s3)
    addarc!(fsm, s3, s4)

    lhs = ones(T, 4, 4)

    label2pdfid = Dict("a" => 1, "b" => 2, "c" => 3, "d" => 4)

    mfsm = MatrixFSM(fsm, label2pdfid)
    mfsm = convert(MatrixFSM{TropicalSemiring{T}}, mfsm)
    μ = maxstateposteriors(mfsm, lhs)
    path = bestpath(mfsm, μ)

    @test join(mfsm.labels[path], " ") == "a b c d"
end

max_stars_repo_name: BUTSpeechFIT/HiddenMarkovModel

max_issues_repo_name: BUTSpeechFIT/HiddenMarkovModel

max_forks_repo_name: BUTSpeechFIT/ProbabilisticFSMs


text:
for N in [Float64, Float32, Rational{Int}]
    H1 = BallInf(N[0, 0], N(1))
    H2 = BallInf(N[-1//10, 0], N(1))
    H3 = Hyperrectangle(N[-16//10, -2], N[1//2, 1//2])

    for p in N[1, 2, Inf]
        @test distance(H1, H2; p=p) == distance(H2, H1; p=p) == N(0)
        @test distance(H1, H3; p=p) == distance(H3, H1; p=p) ≈ norm(N[1//10, 1//2], p)
        @test distance(H2, H3; p=p) == distance(H3, H2; p=p) ≈ norm(N[0, 1//2], p)
    end
end


text:
import CLBLAS
import OpenCL

    const clblas = CLBLAS
    clblas.setup()
    const cl = OpenCL
    alpha = cl.cl_float(5)
    beta = cl.cl_float(2)

    device, context, queue = clblas.get_next_compute_context()
    A = clblas.randBuf(8, 12; ctx=context, queue=queue)
    B = clblas.randBuf(12, 24; ctx=context, queue=queue)
    C = clblas.randBuf(8, 24; ctx=context, queue=queue)

    matA = clblas.fetch(A)
    matB = clblas.fetch(B)
    matC = clblas.fetch(C)

    data = (5.0*matA*matB) + (2.0*matC)
    future = clblas.clblasSgemm(UInt32(0), UInt32(0), alpha, A, B, beta, C)

    result = fetch(future)
    println(result)
    println("------------")
    println(data)
    for i in 1:length(data)
      Test.@test_approx_eq_eps result[i] data[i]  0.00001
    end

    info("success!")
    clblas.teardown()


text:
using ColorSchemes, CairoMakie, FileIO, JuMP
using ProgressMeter

include(joinpath("src", "EMP-model", "model.jl"))
import .GlnModel

imgpath = joinpath("docs", "imgs", "EMP-model", "Concentration_changes")

#=
Warm up problem to get all the symbols in the model.
Create arrays using these symbols as the names.
This is the fast way of doing:
mu = Float64[]
atp = Float64[]
etc.
=#
test_model = GlnModel.gln_model();
syms = Symbol.(all_variables(test_model))
for sym in syms
    eval(:($sym=Float64[]))
end

prev_sol=Dict{Symbol, Float64}() # a variable to store the previous solution

x = []

con_lb = 1e-5
con_ub = 1e-3
@showprogress for var in range(log(con_lb); stop=log(con_ub), length=20)
    lac_ext = exp(var)

    gln_model = GlnModel.gln_model(
        glc_ext = 50e-3,
        lac_ext = lac_ext,
        nh4_ext = 10e-3,
        ac_ext = 1e-16,
        etoh_ext = 1e-6,
        # atp_adp_ratio = -1,
        # nadh_nad_ratio = -1, # unconstrained
        num_ms = 10,
    )

    solve_1 = GlnModel.max_mu!(gln_model, prev_sol)

    if !solve_1
        continue # if mu not found, skip
    end

    solve_2 = GlnModel.find_nearest!(gln_model, prev_sol) # find a solution close to the previous one, smoothen things

    if solve_1 && solve_2
        for sym in syms
            push!(eval(sym), value(gln_model[sym]))
        end
        # update old previous solution
        GlnModel.update!(gln_model, prev_sol)
    end

    push!(x, lac_ext)

end

x = Float64.(x)
y = exp.(atp)./exp.(adp)
z = exp.(nadh)./exp.(nad)

f = Figure()
ax = Axis(f[1,1], xscale=log10,
          xlabel="External Lactate concentration",
          ylabel="Relative ratio change",
)

scatter1 = scatter!(ax, x, y)
line1 = lines!(ax, x, y)
scatter2 = scatter!(ax, x, z)
line2 = lines!(ax, x, z)

f[1, 2] = Legend(
    f,
    [[scatter1, line1], [scatter2, line2]],
    ["ATP/ADP ratio", "NADH/NAD ratio"],
)

f

FileIO.save(joinpath(imgpath, "Ratio_change_Lactate.pdf"), f)


text:
fname = "AliceInWonderland.txt"
fstring = read(fname,String)
farr = split(fstring)

# How many words are there in Alice in Wonderland?
length(farr)

# What is the last word?
farr[end]

# What is the last word before "THE END"?
farr[(length(farr)-10):length(farr)]
farr[length(farr)-2]

# This is an example of how you might count the number of times the word "alice" appears:
count(i->(lowercase(i)=="alice"),farr)

# Write a function that will count the number of times any word appears

function wordcount(x)
    count(i->(lowercase(i)==x),farr)
end

wordcount("the")

# Make the function work even if I passed it an argument with capital letters (e.g. "Alice")?

# Make a collection which tells me how many times the following words appear:
interesting_words = ["alice","mushroom", "caterpillar", "head", "cheshire", "the", "and"]

Wonderland = Dict()
for word in interesting_words
    Wonderland[word] = wordcount(word)
end

print(Wonderland)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Xyce_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Xyce")
JLLWrappers.@generate_main_file("Xyce", UUID("b812464d-768c-53cf-a694-ea0631d6fbf4"))
end  # module Xyce_jll


text:
__precompile__()
"""
    EasyEcon

Easy & friendly mathematics/economics in OLG model;
recommend to using it for fast debugging & mass use
"""
module EasyEcon

# ==============================================================================
    """
        CDProdFunc( A::Real, K::Real, L::Real, β::Real; κ::Real = 0.05, GetPrice::Bool = true, TechType::String = "Hicks"  )

    Cobb-Douglas production function, where A is technology, K is capital, L is labour, β is income share of capital;
    allows three types of technology: "Hicks", "Harod" & "Solow", controled by TechType parameter;
    if GetPrice is true, separately returns: output (Y), equilibrium NET interest rate (r-κ) & equilibrium wage level (w̄);
    where the capital depreciation rate is controled by keyword parameter κ;
    if GetPrice is false, only returns: output (Y), a real number
    """
    function CDProdFunc( A::Real, K::Real, L::Real, β::Real; κ::Real = 0.05, GetPrice::Bool = true, TechType::String = "Hicks"  )
        # 1. bounding check
        tmp = [A,K,L,β,κ]
        @assert( all(tmp .> 0), string("bound failure [A,K,L,β,κ], all positive required: ",tmp) )
        # 2. compute
        if TechType == "Hicks"
            Y = A * K^β * L^(1-β)
            if GetPrice
                r = β * A * (K/L)^(β-1) - κ
                w̄ = (1-β) * A * (K/L)^β
                return Y::Real, r::Real, w̄::Real
            else
                return Y::Real
            end
        elseif TechType == "Harod"
            Y = K^β * (A * L)^(1-β)
            if GetPrice
                r = β * (K/A/L)^(β-1) - κ
                w̄ = (1-β) * A * (K/A/L)^β
                return Y::Real, r::Real, w̄::Real
            else
                return Y::Real
            end
        elseif TechType == "Solow"
            Y = (A * K)^β * L^(1-β)
            if GetPrice
                r = β * (A*K/L)^(β-1) - κ
                w̄ = (1-β) * (A*K/L)^β
                return Y::Real, r::Real, w̄::Real
            else
                return Y::Real
            end
        else
            throw( ErrorException("invalid TechType, pls choose from \"Hicks\", \"Harod\", \"Solow\"") )
        end
        return nothing
    end

    # --------------------------------------------
    """
        WageProfile( w̄::Real, ε::Vector, Nst::Vector, Lst::Vector; GetScalingCoef::Bool = false  )

    profile wage among different working cohorts in year t where age ``s=1,...,S_r`` (``S_r`` is retiring age);
    using equation: ``w̄ L_{t} = \\sum^{S_r}_{s=1} w_{s,t} L_{s,t} N_{s,t}``;
    where ``L_{t} = \\sum^{S_r} L_{s,t}N_{s,t}`` is labour supply, ``w_{s,t} = o_{t} \varepsilon_{s}`` is wage profiling,
    ``o_{t}`` is a year-specific scaling coefficient (computed by the equation), ``\varepsilon_{s}`` is age-specific profiling coefficients;
    using Lst for ``L_{s,t}``, and Nst for ``N_{s,t}``;
    if GetSaclingCoef is true, separately returns: wprofiled::Vector (``w_{s,t}``), o::Real (``o_{t}``);
    if false, only returns: wprofiled;

    ## Inputs:
    1. ̄w: average wage level
    2. ε: wage profiling coefficient
    3. Nst: working population vector
    4. Lst: labour vector (1 - leisure)
    5. GetScalingCoef: if to return the wage scaling coefficient
    """
    function WageProfile( w̄::Real, ε::Vector, Nst::Vector, Lst::Vector; GetScalingCoef::Bool = false  )
        # 1. size check
        @assert( length(ε) == length(Nst) == length(Lst), "requires same length: ε, Nst, Lst" )
        # 2. compute
        L = sum( Nst .* Lst )  # get total labour supply
        o = L / sum( Nst .* Lst .* ε )  # get scaling coefficient
        wprofiled = w̄ .* o .* ε  # get profiled wage levels
        # 3. return
        if GetScalingCoef
            return wprofiled::Vector, o::Real
        else
            return wprofiled::Vector
        end
        return nothing
    end

    # --------------------------------------------
    """
        PAYGPension( πCoef::Real, w::Vector, N::Vector, Lab::Vector, Sr::Int )

    Computes an average pension benefit amount for retired popultion:
        `` \\pi_{t} \\sum^{S_r}_{s=1} w_{s,t} N_{s,t} L_{s,t} = \\Lambda_{t} \\sum^{S}_{s>S_r} N_{s,t}``
    where ``\\pi_{t}`` is total contribution rate in year ``t``, ``L_{s,t}`` is labor supply vector, ``\\Lambda_{t}`` is average pension benefit amount in year ``t``, ``N_{s,t}`` is population.

    ##Inputs:
    1. πCoef : ``\\pi_{t}``
    2. w : ``w_{s,t}``, nomial wage level, len = Sr
    3. N : compelete population, from 1 to Sr (retirement age) to S (max age), len = S
    4. Lab : ``L_{s,t} = 1 - l_{s,t}``, where ``l_{s,t}`` is leisure, len = Sr
    5. Sr : retirement age

    returns a Real number, the average pension benefit amount ``\\Lambda_{t}``.
    """
    function PAYGPension( πCoef::Real, w::Vector, N::Vector, Lab::Vector, Sr::Int )
        # 1. size check
        @assert( length(w) == length(N[1:Sr]) == length(Lab), "uncompatible vector size: w, Nwork or Lab"  )
        # 2. accounting
        local Λ::Float64
        Λ = πCoef * sum( w .* N[1:Sr] .* Lab ) / sum( N[Sr+1:end] )
        return Λ
    end

    # --------------------------------------------
    """
        Get𝕡( 𝕓::Real, πMf::Real, w::Vector, L::Vector, N::Vector, Sr::Int )

    Gets transfer payment amount 𝕡 from firm contribution to UE-BMI in year ``t`` to those retired in this year.
    based on a transfer rate 𝕓.

    ## Inputs:
    1. 𝕓: transfer rate/fraction from firm total contribution amount (to UE-BMI) to retired generations
    2. πMf: total contribution rate by firm on nomial wage level ``w_{s,t}``; it is defined as: ``\\pi^{Mf} = \\frac{\\zeta_{t}}{1+z \\eta + \\zeta}``
    3. w: nomial wage level vector, len = Sr
    4. Lab: labor supply, len = Sr
    5. N: population, len = S (env.MAX_AGE)
    6. Sr: retirement age

    returns a Real number 𝕡, the transfer payment amount per capita to retired generations.
    """
    function Get𝕡( 𝕓::Real, πMf::Real, w::Vector, Lab::Vector, N::Vector, Sr::Int )
        # 1. size check
        @assert( length(w) == length(Lab) == length(N[1:Sr]) , "uncompatible vector size: w, Lab or N" )
        # 2. accounting
        local 𝕡::Float64
        𝕡 = 𝕓 * πMf * sum( w .* N[1:Sr] .* Lab ) / sum( N[Sr+1:end] )
        return 𝕡
    end

# ==============================================================================
end  # module ends


text:
using Documenter, GalacticOptim

makedocs(
    sitename="GalacticOptim.jl",
    authors="Chris Rackauckas, Vaibhav Kumar Dixit et al.",
    clean=true,
    doctest=false,
    modules=[GalacticOptim],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://galacticoptim.sciml.ai/stable/"),

    pages=[
        "GalacticOptim.jl: Unified Global Optimization Package" => "index.md",

        "Tutorials" => [
            "Basic usage" => "tutorials/intro.md",
            "Rosenbrock function" => "tutorials/rosenbrock.md",
            "Minibatch" => "tutorials/minibatch.md",
            "Symbolic Modeling" => "tutorials/symbolic.md"
        ],

        "API" => [
            "OptimizationProblem" => "API/optimization_problem.md",
            "OptimizationFunction" => "API/optimization_function.md",
            "solve" => "API/solve.md",
            "ModelingToolkit Integration" => "API/modelingtoolkit.md"
        ],
        "Optimizer Packages" => [
            "BlackBoxOptim.jl" => "optimization_packages/blackboxoptim.md",
            "CMAEvolutionStrategy.jl" => "optimization_packages/cmaevolutionstrategy.md",
            "Evolutionary.jl" => "optimization_packages/evolutionary.md",
            "Flux.jl" => "optimization_packages/flux.md",
            "GCMAES.jl" => "optimization_packages/gcmaes.md",
            "MathOptInterface.jl" => "optimization_packages/mathoptinterface.md",
            "MultistartOptimization.jl" => "optimization_packages/multistartoptimization.md",
            "Metaheuristics.jl" => "optimization_packages/metaheuristics.md",
            "NOMAD.jl" => "optimization_packages/nomad.md",
            "NLopt.jl" => "optimization_packages/nlopt.md",
            "Nonconvex.jl" => "optimization_packages/nonconvex.md",
            "Optim.jl" => "optimization_packages/optim.md",
            "QuadDIRECT.jl" => "optimization_packages/quaddirect.md"
        ],
    ]
)

deploydocs(
    repo="github.com/SciML/GalacticOptim.jl";
    push_preview=true
)


text:
####
#### OpSum
####

struct OpSum{OpT<:AbstractOp, StringT, CoeffT} <: AbstractSum{OpT, StringT, CoeffT}
    strings::StringT
    coeffs::CoeffT

    function OpSum(strings, coeffs; already_sorted=false)
        if length(strings) != length(coeffs)
            throw(DimensionMismatch("bad dims"))
        end
        if ! isempty(strings)
            if ! already_sorted  # Slightly dangerous to only do length check if not sorted
                n = length(first(strings))
                if ! all(x -> length(x) == n, strings)
                    throw(DimensionMismatch("Operator strings are of differing lengths."))
                end
                sort_and_sum_duplicates!(strings, coeffs)
            end
        end
        return new{eltype(eltype(strings)), typeof(strings), typeof(coeffs)}(strings, coeffs)
    end
end

strip_typeof(::OpSum{W, T, CoeffT}) where {W, T, CoeffT} = OpSum{W}

"""
    term_type(::Type{<:OpSum{T}})

Return the the concrete `OpTerm` type associated with `OpSum{T}`.
This is used in generic code to construct terms given the type of a sum.
"""
term_type(::Type{<:OpSum{T}}) where T = OpTerm{T}

"""
    sum_type(::Type{<:OpTerm{T}})

Return the the concrete `OpSum` type associated with `OpTerm{T}`.
This is used in generic code to construct sums given the type of a term.
"""
sum_type(::Type{<:OpTerm{T}}) where T = OpSum{T}

#OpSum{T}(strings, coeffs; already_sorted=false) where T = OpSum(strings, coeffs; already_sorted=already_sorted)
#OpSum{T}(args...) where T = OpSum(args...)
OpSum{T}(args...; kwargs...) where T = OpSum(args...; kwargs...)

function OpSum{T}(strings::AbstractVector{<:AbstractString}, coeffs) where T
    return OpSum(Vector{T}.(strings), coeffs)
end

function OpSum{T}(strings::AbstractVector{<:AbstractString}) where T
    coeffs = fill(_DEFAULT_COEFF, length(strings))
    return OpSum{T}(strings, coeffs)
end

function OpSum(v::AbstractVector{<:OpTerm}; already_sorted=false)
    strings = [op_string(x) for x in v]
    coeffs = [x.coeff for x in v]
    return OpSum(strings, coeffs; already_sorted=already_sorted)
end

"""
    empty(s::OpSum)

Create an empty `OpSum`, that is with no terms, similar to `s`.
"""
Base.empty(s::OpSum) = strip_typeof(s)(empty(s.strings), empty(s.coeffs))

Base.isempty(s::OpSum) = isempty(s.strings) && isempty(s.coeffs)  # Do we need both. Should we check ?

"""
    OpSum(v::AbstractMatrix{<:AbstractOp}, coeffs=fill(_DEFAULT_COEFF, size(v, 1)))

Construct a sum from a matrix of single-particle operators. If `size(v) == (m, n)`, then
the the sum has `m` terms with `n` operators in each string.
"""
function OpSum(v::AbstractMatrix{<:AbstractOp}, coeffs=fill(_DEFAULT_COEFF, size(v, 1)))
    strings = @inbounds [v[i,:] for i in 1:size(v,1)]
    return OpSum(strings, coeffs)
end

OpSum{OpT}() where OpT <: AbstractOp = OpSum(Vector{OpT}[], Complex{Float64}[])

"""
    rand_op_sum(::Type{OpT}, n_factors::Integer, n_terms::Integer; coeff_func=nothing) where {OpT <: AbstractOp}

Return an `OpSum{OpT}` of `n_terms` terms of `n_factors` factors each.

If `coeff_func` is `nothing`, then the coefficients are all equal to one. Otherwise `coeff_func` must be a function that
takes one argument `n_terms`, and returns `n_terms` coefficients.

# Examples
```julia
julia> rand_op_sum(Pauli, 4, 3)
(1 + 0im) * IIYX
(1 + 0im) * IIZY
(1 + 0im) * ZXIX

julia> rand_op_sum(FermiOp, 3, 4; coeff_func=randn)
N-I * -1.125005096910286
EI- * 0.8400663971914751
EE- * -0.11901005957119838
-E- * 0.4507511196805674
```
"""
function rand_op_sum(::Type{OpT}, n_factors::Integer, n_terms::Integer; coeff_func=nothing) where {OpT <: AbstractOp}
    term_string = [rand(OpT, n_factors) for i in 1:n_terms]
    if coeff_func === nothing
        coeffs = fill(_DEFAULT_COEFF, n_terms)
    else
        coeffs = coeff_func(n_terms)
    end
    return OpSum{OpT}(term_string, coeffs)
end

#####
##### Conversion
#####

"""
    Matrix(ps::OpSum)

Convert `ps` to a dense `Matrix`.

# Examples

We convert a matrix to a `OpSum` and then back to a matrix.
```jldoctest
julia> m = [0.1 0.2; 0.3 0.4];

julia> PauliSum(m)
(0.25 + 0.0im) * I
(0.25 + 0.0im) * X
(0.0 - 0.04999999999999999im) * Y
(-0.15000000000000002 + 0.0im) * Z

julia> Matrix(PauliSum(m)) ≈ m
true
```
"""
Base.Matrix(ps::OpSum) = Matrix(SparseArrays.sparse(ps))

## ThreadsX helps enormously for large sums. 22x faster for 4^8 terms
## But, ThreadsX is five times slower for 1-q sum
"""
    sparse(ps::OpSum)

Convert `ps` to a sparse matrix.
"""
function SparseArrays.sparse(ps::OpSum)
    if length(ps) >= 64
        return ThreadsX.sum(SparseArrays.sparse(ps[i]) for i in eachindex(ps))
    else
        return sum(SparseArrays.sparse(ps[i]) for i in eachindex(ps))
    end
end

####
#### Compare / predicates
####

## TODO: improve this. Check for coefficients, etc.
function Base.iszero(opsum::OpSum)
    return iszero(length(opsum)) || all(iszero, opsum.coeffs)
end

Base.zero(opsum::OpSum) = empty(opsum)

function ZChop.zchop(opsum::OpSum, zeps::Real=ZChop.ZEPS)
    newopsum = empty(opsum)
    newstrings = newopsum.strings
    newcoeffs = newopsum.coeffs
    for i in eachindex(opsum.coeffs)
        c = opsum.coeffs[i]
        if ! iszero(ZChop.zchop(c, zeps))
            push!(newstrings, opsum.strings[i])
            push!(newcoeffs, c)
        end
    end
    return newopsum
end

function ZChop.zchop!(opsum::OpSum, zeps::Real=ZChop.ZEPS)
    inds = Int[]
    for i in eachindex(opsum.coeffs)
        if iszero(ZChop.zchop(opsum.coeffs[i], zeps))
            push!(inds, i)
        end
    end
    deleteat!(opsum.coeffs, inds)
    deleteat!(opsum.strings, inds)
    return opsum
end

####
#### Algebra / mathematical operations
####

# We use lmul! because that's how LinearAlgebra offers "scaling" of a Matrix (or rmul!)
"""
    lmul!(psum::OpSum, n)

Left multiplies the coefficient of `psum` by `n` in place.
"""
function LinearAlgebra.lmul!(opsum::OpSum, n)
    coeffs = opsum.coeffs
    # FastBroadcast.@.. opsum.coeffs = n * opsum.coeffs
    # @. opsum.coeffs = n * opsum.coeffs
    @inbounds for i in eachindex(coeffs)
        coeffs[i] = coeffs[i] * n
    end
    return opsum
end


text:
module NLSolvers

using Gridap
using NLsolve

export NLSolver

import Gridap: solve!

mutable struct NLSolver <: NonLinearSolver
  ls::LinearSolver
  kwargs::Dict
end

function NLSolver(ls::LinearSolver;kwargs...)
  @assert ! haskey(kwargs,:linsolve) "linsolve cannot be used here. It is managed internally"
  NLSolver(ls,kwargs)
end

function NLSolver(;kwargs...)
  ls = BackslashSolver()
  NLSolver(ls;kwargs...)
end

export JuliaNLSolver

@deprecate JuliaNLSolver(;kwargs...) NLSolver(;kwargs...)

@deprecate JuliaNLSolver(ls;kwargs...) NLSolver(ls;kwargs...)

mutable struct NLSolversCache
  df::OnceDifferentiable
  ss::SymbolicSetup
  result
end

function solve!(x::AbstractVector,nls::NLSolver,op::NonLinearOperator)
  cache = _setup_cache(x,nls,op)
  solve!(x,nls,op,cache)
  cache
end

function solve!(
  x::AbstractVector,nls::NLSolver,op::NonLinearOperator,cache::NLSolversCache)
  df = cache.df
  ss = cache.ss
  kwargs = nls.kwargs
  function linsolve!(x,A,b)
    ns = numerical_setup(ss,A)
    solve!(x,ns,b)
  end
  r = nlsolve(df,x;linsolve=linsolve!,kwargs...)
  cache.result = r
  x[:] .= r.zero
end

function _setup_cache(x0,nls,op)

  f!(r,x) = residual!(r,op,x)
  j!(j,x) = jacobian!(j,op,x)

  f0 = residual(op,x0)
  j0 = jacobian(op,x0)

  df = OnceDifferentiable(f!,j!,x0,f0,j0)

  ss = symbolic_setup(nls.ls,j0)

  NLSolversCache(df,ss,nothing)
end

import Base: show

function show(io::IO,self::NLSolversCache)
  print(io,"$(nameof(typeof(self))) object")
end

end # module


text:
#!/usr/bin/env julia
using InvariantCausal
using Test
println("Starting tests")
@time begin
include("test_regression.jl")
include("test_search.jl")
end

max_issues_repo_name: richardkwo/InvariantCausal

max_forks_repo_name: richardkwo/InvariantCausal


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using DelimitedFiles

mmesh = Mesh3(10, 10, 10, 1/10)

u = readdlm("val.txt")[:]
visualize_scalar_on_fem_points(u, mmesh, filename = "laplace.png")


text:
using OrdinaryDiffEq, ForwardDiff
import OrdinaryDiffEq.ODEProblem

export ContinuousDS, ODEProblem

#######################################################################################
#                                     Constructors                                    #
#######################################################################################
"""
    ContinuousDS(state, eom [, jacob]) <: DynamicalSystem
`D`-dimensional continuous dynamical system (used for `D ≤ 10`).
# Fields:
* `state::SVector{D}` : Current state-vector of the system, stored in the data format
  of `StaticArray`'s `SVector`.
* `eom::F` (function) : The function that represents the system's equations of motion
  (also called vector field). The function is of the format: `eom(u) -> SVector`
  which means that given a state-vector `u` it returns an `SVector` containing the
  derivatives `du` at the current state.
* `jacob::J` (function) : A function that calculates the system's jacobian matrix,
  based on the format: `jacob(u) -> SMatrix` which means that given a state-vector
  `u` it returns an `SMatrix` containing the Jacobian at that state.
  If the `jacob` is not provided by the user, it is created with *tremendous* efficiency
  using the module `ForwardDiff`. Most of the time, for low dimensional systems, this
  Jacobian is within a few % of speed of a user-defined one.
"""
mutable struct ContinuousDS{D, T<:Real, F, J} <: DynamicalSystem
  state::SVector{D,T}
  eom::F
  jacob::J
end
# constsructor without jacobian (uses ForwardDiff)
function ContinuousDS(u0::AbstractVector, eom)
  su0 = SVector{length(u0)}(u0)
  @inline ForwardDiff_jac(x) = ForwardDiff.jacobian(eom, x)
  return ContinuousDS(su0, eom, ForwardDiff_jac)
end
function ContinuousDS(u0::AbstractVector, eom, jac)
  su0 = SVector{length(u0)}(u0)
  return ContinuousDS(su0, eom, jac)
end

dimension(::ContinuousDS{D, T, F, J})  where {D<:ANY, T<:ANY, F<:ANY, J<:ANY} = D
#######################################################################################
#                                 System Evolution                                    #
#######################################################################################
"""
```julia
ODEProblem(ds::ContinuousDS, t)
```
Return a type `ODEProblem` with the given
system information (t0 is zero).
This can be passed directly into `solve` from `DifferentialEquations`.
"""
function ODEProblem(ds::ContinuousDS, t)
  odef = (t, u) -> ds.eom(u)
  OrdinaryDiffEq.ODEProblem(odef, ds.state, (zero(t), t))
end

function get_sol(prob::ODEProblem, diff_eq_kwargs::Dict = Dict())
  if haskey(diff_eq_kwargs, :solver)
    solver = diff_eq_kwargs[:solver]
    pop!(diff_eq_kwargs, :solver)
    sol = solve(prob, solver; diff_eq_kwargs..., save_everystep=false)
  else
    sol = solve(prob, Tsit5(); diff_eq_kwargs..., save_everystep=false)
  end
  return sol
end

function evolve(ds::ContinuousDS, t::Real = 1.0; diff_eq_kwargs::Dict=Dict())
  prob = ODEProblem(ds, t)
  return get_sol(prob, diff_eq_kwargs)[end]
end
function evolve(
  state::AbstractVector, ds::ContinuousDS, t::Real=1.0;
  diff_eq_kwargs::Dict=Dict())

  prob = ODEProblem(ds, t)
  prob.u0 = state
  return get_sol(prob, diff_eq_kwargs)[end]
end

function evolve!(ds::ContinuousDS, t::Real = 1.0; diff_eq_kwargs::Dict = Dict())
  ds.state = evolve(ds, t, diff_eq_kwargs = diff_eq_kwargs)
  return ds
end

"""
```julia
evolve!(prob::ODEProblem [, t::Real, diff_eq_kwargs::Dict]) -> final_state
```
Evolve the problem using the solvers of `DifferentialEquations`, update
the problem's state as the final state of the solution and return that state.

If `t` is given, the problem is evolved for that much time (else the existing
`tspan` is used). Notice that in this function, `diff_eq_kwargs` is *not* a keyword
argument.
"""
function evolve!(prob::ODEProblem, t::Real, diff_eq_kwargs::Dict)
  prob.tspan = (zero(t), t)
  if haskey(diff_eq_kwargs, :solver)
    solver = diff_eq_kwargs[:solver]
    pop!(diff_eq_kwargs, :solver)
    state = solve(prob, solver; diff_eq_kwargs..., save_everystep=false)[end]
  else
    state = solve(prob, Tsit5(); diff_eq_kwargs..., save_everystep=false)[end]
  end
  prob.u0 = state
  return state
end
function evolve!(prob::ODEProblem, diff_eq_kwargs::Dict)
  if haskey(diff_eq_kwargs, :solver)
    solver = diff_eq_kwargs[:solver]
    pop!(diff_eq_kwargs, :solver)
    state = solve(prob, solver; diff_eq_kwargs..., save_everystep=false)[end]
  else
    state = solve(prob, Tsit5(); diff_eq_kwargs..., save_everystep=false)[end]
  end
  prob.u0 = state
  return state
end
function evolve!(prob::ODEProblem)
  state = solve(prob, Tsit5(); save_everystep=false)[end]
  prob.u0 = state
  return state
end

function timeseries(ds::ContinuousDS, T::Real;
  dt::Real=0.05, diff_eq_kwargs::Dict=Dict(), mutate = true)

  T<=0 && throw(ArgumentError("Total time `T` must be positive."))
  D = dimension(ds)
  t = zero(T):dt:T
  prob = ODEProblem(ds, T)
  kw = Dict{Symbol, Any}(diff_eq_kwargs)
  kw[:saveat] = t
  sol = get_sol(prob, kw)
  TS = zeros(eltype(ds.state), length(t), D)
  for j in 1:D
    TS[:, j] .= sol[j,:]
  end
  # using: transpose(hcat(get_sol(prob, diff_eq_kwargs).u...))
  # is so absurdly tragically slower.
  if mutate
    ds.state = TS[end, :]
  end
  return TS
end

#######################################################################################
#                                 Tangent Space                                       #
#######################################################################################
function tangentbundle_setup(ds::ContinuousDS, dt)
  D = dimension(ds)
  S = [ds.state eye(eltype(ds.state), D)]
  f = ds.eom
  jac = ds.jacob
  function tbeom(t, u, du)
    du[:, 1] .= f(u)
    A_mul_B!(view(du, :, 2:D+1), jac(view(u, :, 1)), view(u, :, 2:D+1))
  end
  tbprob = ODEProblem(tbeom, S, (zero(dt), dt))
  return tbprob
end

#######################################################################################
#                                 Pretty-Printing                                     #
#######################################################################################
import Base.show
function Base.show(io::IO, s::ContinuousDS{N, S, F, J}) where {N<:ANY, S<:ANY, F<:ANY, J<:ANY}
  print(io, "$N-dimensional continuous dynamical system:\n",
  "state: $(s.state)\n", "e.o.m.: $F\n", "jacobian: $J")
end

@require Juno begin
  function Juno.render(i::Juno.Inline, s::ContinuousDS{N, S, F, J}) where {N<:ANY, S<:ANY, F<:ANY, J<:ANY}
    t = Juno.render(i, Juno.defaultrepr(s))
    t[:head] = Juno.render(i, Text("$N-dimensional continuous dynamical system"))
    t
  end
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Deep learning with Flux")

#

# # Recognizing handwritten digits using a neural network

#-

# We have now reached the point where we can tackle a very interesting task: applying the knowledge we have gained with machine learning in general, and `Flux.jl` in particular, to create a neural network that can recognize handwritten digits! The data are from a data set called MNIST, which has become a classic in the machine learning world.
#
# [We could also try to apply the techniques to the original images of fruit instead. However, the fruit images are much larger than the MNIST images, which makes the learning a suitable neural network too slow.]

#-

# ## Data munging

#-

# As we know, the first difficulty with any new data set is locating it, understanding what format it is stored in, reading it in and decoding it into a useful data structure in Julia.
#
# The original MNIST data is available [here](http://yann.lecun.com/exdb/mnist); see also the [Wikipedia page](https://en.wikipedia.org/wiki/MNIST_database). However, the format that the data is stored in is rather obscure.
#
# Fortunately, various packages in Julia provide nicer interfaces to access it. We will use the one provided by `Flux.jl`.
#
# The data are images of handwritten digits, and the corresponding labels that were determined by hand (i.e. by humans). Our job will be to get the computer to **learn** to recognize digits by learning, as usual, the function that relates the input and output data.

#-

# ### Loading and examining the data

#-

# First we load the required packages:

using Flux, Flux.Data.MNIST, Images

# Now we read in the data:

labels = MNIST.labels();
images = MNIST.images();

# We see that `images` is a `Vector`, i.e. an `Array{T, 1}` with a complicated parameter `T`. It has length

length(labels)

# But we can just look at the first handful to get a sense of the contents

images[1:5]

#-

labels[1:5]' # transposed to match the above

# So the $i$th entry of the array is the data for the $i$th image.

typeof(images[1])

# As with the fruit images from the start of the course, the image is an array of color blocks, except that now each pixel just has a grey scale.
#
# To see the actual numeric content of the image, we can do, for example

Float64.(images[1])

# ### Setting up a neural network

#-

# In the previous notebooks, we arranged the input data for Flux as a `Vector` of `Vector`s.
# Now we will use an alternative arrangement, as a matrix, since that allows `Flux` to use matrix operations, which are more efficient.
#
# The column $i$ of the matrix is a vector consisting of the $i$th data point $\mathbf{x}^{(i)}$.  Similarly, the desired outputs are given as a matrix, with the $i$th column being the desired output $\mathbf{y}^{(i)}$.

n_inputs = unique(length.(images))[]

#-

n_outputs = length(unique(labels))

# #### Creating the features
#
# We want to create a vector of feature vectors, each with the floating point values of the 784 pixels.
#
# An image is a matrix of colours, but now we need a vector of floating point numbers instead. To do so, we just arrange all of the elements of the matrix in a certain way into a single list; fortunately, Julia already provides the function `vec` to do so!

#-

# Let's use a subset of $N=5,000$ of the total $60,000$ images available in order to hold out test images that our model hasn't been trained on.

preprocess(img) = vec(Float64.(img))

#-

xs = preprocess.(images[1:5000]);

# #### Creating the labels
#
# We can just use `Flux.onehot` to create them:

ys = [Flux.onehot(label, 0:9) for label in labels[1:5000]];

# #### Create the batched matrices for efficiency

#-

# Let's also create a function so we can easily create independent batches from arbitrary segments of the original dataset.

function create_batch(r)
    xs = [preprocess(img) for img in images[r]]
    ys = [Flux.onehot(label, 0:9) for label in labels[r]]
    return (Flux.batch(xs), Flux.batch(ys))
end

# We'll train our model on the first 5000 images.

trainbatch = create_batch(1:5000);

# ## Setting up the neural network

#-

# Now we must set up a neural network. Since the data is complicated, we may expect to need several layers.
# But we can start with a single layer.
#
# - The network will take as inputs the vectors $\mathbf{x}^{(i)}$, so the input layer has $n$ nodes.
#
# - The output will be a one-hot vector encoding the digit from 1 to 9 or 0 that is desired. There are 10 possible categories, so we need an output layer of size 10.
#
# It is then our task as neural network designers to insert layers between these input and output layers, whose weights will be tuned during the learning process. *This is an art, not a science*! But major advances have come from finding a good structure for the network.

model = Chain(Dense(n_inputs, n_outputs, identity), softmax)
L(x,y) = Flux.crossentropy(model(x), y)
opt = SGD(params(model))
@time Flux.train!(L, [trainbatch], opt)
@time Flux.train!(L, [trainbatch], opt)

# ## Training

#-

# Just as before, use `repeated` to create an **iterator**. It does not copy the data 100 times, which would be very wasteful; it just gives an object that repeatedly loops over the same data:

Iterators.repeated(trainbatch, 100);

# We can see what the total current loss is (after training just a handful to times above):

L(trainbatch...)

#-

@time Flux.train!(L, [trainbatch], opt)

#-

L(trainbatch...)

# ### Using callbacks

#-

# The `train!` function can take an optional keyword argument, `cb` (short for "**c**all**b**ack"). A callback function is a function that you provide as an argument to a function `f`, which "calls back" your function every so often.
#
# This provides the possibility to provide a function that is called at each step or every so often during the training process.
# A common use case is to provide a visual trace of the training process by printing out the current value of the `loss` function:

callback() = @show(L(trainbatch...))

Flux.train!(L, Iterators.repeated(trainbatch, 3), opt; cb = callback)

# However, it is expensive to calculate the complete `loss` function and it is not necessary to output it every step. So `Flux` also provides a function `throttle`, that provides a mechanism to call a given function at most once every certain number of seconds:

Flux.train!(L, Iterators.repeated(trainbatch, 40), opt; cb = Flux.throttle(callback, 1))

# Of course, that's just measuring the loss over the same data it's training on. It'd be more representative to test against novel data. In fact, let's track the performance of both as we continue training our model. In order to do so, we need to create a batch of test data.

testbatch = create_batch(5001:10000);

#-

using Printf
train_loss = Float64[]
test_loss = Float64[]
function update_loss!()
    push!(train_loss, L(trainbatch...).data)
    push!(test_loss, L(testbatch...).data)
    @printf("train loss = %.2f, test loss = %.2f\n", train_loss[end], test_loss[end])
end

#-

Flux.train!(L, Iterators.repeated(trainbatch, 1000), opt;
                cb = Flux.throttle(update_loss!, 1))

#-

using Plots
plot(1:length(train_loss), train_loss, xlabel="~seconds of training", ylabel="loss", label="train")
plot!(1:length(test_loss), test_loss, label="test")

# ## Testing phase

#-

# We now have trained a model, i.e. we have found the parameters `W` and `b` for the network layer(s). In order to **test** if the learning procedure was really successful, we check how well the resulting trained network performs when we test it with images that the network has not yet seen!
#
# Often, a dataset is split up into "training data" and "test (or validation) data" for this purpose, and indeed the MNIST data set has a separate pool of training data. We can instead use the images that we have not included in our reduced training process.

i = 5001
display(images[i])
labels[i], findmax(model(preprocess(images[i]))) .- (0, 1)

#-

model(preprocess(images[i]))

# ## Evaluation
#
# What percent of images are we correctly classifying if we take the highest element to be the chosen answer?

prediction(i) = findmax(model(preprocess(images[i])))[2]-1 # returns (max_value, index)

#-

sum(prediction(i) == labels[i] for i in 1:5000)/5000

#-

sum(prediction(i) == labels[i] for i in 5001:10000)/5000

# ## Improving the prediction

#-

# So far we have used a single layer. In order to improve the prediction, we probably need to use more layers. Try adding more layers yourself and see how the performance changes.

n_hidden = 20
model = Chain(Dense(n_inputs, n_hidden, relu),
              Dense(n_hidden, n_outputs, identity), softmax)
L(x,y) = Flux.crossentropy(model(x), y)
opt = ADAM(params(model))

#-

train_loss = Float64[]
test_loss = Float64[]
Flux.train!(L, Iterators.repeated(trainbatch, 1000), opt;
            cb = Flux.throttle(update_loss!, 1))

#-

plot(1:length(train_loss), train_loss, xlabel="~seconds of training", ylabel="loss", label="train")
plot!(1:length(test_loss), test_loss, label="test")

# ## What about image structure?

#-

# As a final note, notice that our model doesn't take into account any aspect of the image's connected-ness.

using Random
p = randperm(28)
images[1][p,p]

#-

show(preprocess(images[1]))

#-

#nb ?Conv
#jl @doc Conv

max_stars_repo_path: Courses/Deep learning with Flux/1900.Recognizing-handwriting-using-a-neural-network.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
immutable GradientFunction{T <: Function}
    f::T
    name::String
end

grad_funcs = (
    GradientFunction(
        FiniteDiff.gradient,
        "FiniteDiff.gradient(f, x)",
    ),
    GradientFunction(
        (f, x) -> FiniteDiff.gradient(f, x, FiniteDiff.ForwardMode()),
        "FiniteDiff.gradient(f, x, FiniteDiff.ForwardMode())",
    ),
    GradientFunction(
        (f, x) -> FiniteDiff.gradient(f, x, FiniteDiff.BackwardMode()),
        "FiniteDiff.gradient(f, x, FiniteDiff.BackwardMode())",
    ),
    GradientFunction(
        (f, x) -> FiniteDiff.gradient(f, x, FiniteDiff.CentralMode()),
        "FiniteDiff.gradient(f, x, FiniteDiff.CentralMode())",
    ),
    GradientFunction(
        (f, x) -> FiniteDiff.gradient(f, x, FiniteDiff.ComplexMode()),
        "FiniteDiff.gradient(f, x, FiniteDiff.ComplexMode())",
    ),
    GradientFunction(
        Calculus.finite_difference,
        "Calculus.finite_difference(f, x)",
    ),
)


text:
using PowerModelsDistribution
using PowerModels
using Ipopt
using JuMP
using DataFrames
using CSV

PowerModels.logger_config!("error")

include("post_flex_model_new.jl")

data_filename = "data/123Bus_GMLC/IEEE123Master.dss"
data = PowerModelsDistribution.parse_file(data_filename)
filename = "data/ADC_Placement_by_Voltage_Drop.csv"
# filename_flex = "data/ieee123_flex.csv"

filename_flex = ARGS[1]
# data_new = data_surgery(data,filename)
data_new = data_surgery(data, filename, filename_flex)
# result = run_ac_tp_opf_lm(data, with_optimizer(Ipopt.Optimizer))

# model = PowerModels.build_model(data,ACPPowerModel,post_tp_opf_lm_mod)
# result = PowerModels.optimize_model!(model,with_optimizer(Ipopt.Optimizer))
# result_old = _PMs.run_model(data, ACPPowerModel, with_optimizer(Ipopt.Optimizer), post_tp_opf_lm_mod; multiconductor=true, ref_extensions=[ref_add_arcs_trans!])

pm_new = _PMs.build_model(data_new, ACPPowerModel, post_tp_opf_lm_mod; multiconductor=true, ref_extensions=[ref_add_arcs_trans!])
result_new = _PMs.optimize_model!(pm_new, with_optimizer(Ipopt.Optimizer))

ADC = build_ADC_list(filename)
ADC_pm = convert_ADC_pm(data,ADC)
pm_to_source, source_to_pm = build_bus_renaming_dict(data_new)
bus_loads = build_bus_loads(data_new)

to_write = Dict{String,Array{Float64,1}}()
for adc in keys(ADC)
    if !(adc in keys(source_to_pm))
        continue
    end
    # @show adc
    adc_pm = source_to_pm[adc]
    load = bus_loads[adc_pm][1]
    load_id = parse(Int64,load)
    p = zeros(3)
    q = zeros(3)
    for i=1:3
        p[i] = JuMP.value(PowerModels.var(pm_new, 0, i, :pd)[load_id])*1e5
        q[i] = JuMP.value(PowerModels.var(pm_new, 0, i, :qd)[load_id])*1e5
    end
    to_write[adc] = [sum(p),sum(q)]
end

df = DataFrame(ADC = Int[], P = Float64[], Q = Float64[])
for (adc, setpoint) in to_write
    push!(df, Dict(:ADC => parse(Int64,adc), :P => setpoint[1], :Q => setpoint[2]) ) # [parse(Int64,adc), parse)Float64,setpoint[1]), parse(Float64,setpoint[2])])
end

CSV.write("PFO_output.csv", df)

max_stars_repo_name: areiman/PFO-ADC-DER-Testbed

max_issues_repo_name: areiman/PFO-ADC-DER-Testbed

max_forks_repo_name: areiman/PFO-ADC-DER-Testbed


text:
# # Moving sphere in an infinite fluid

# Source code: [`sphere_dipole_tut.jl`](sphere_dipole_tut.jl)

# ## Description

# A rigid sphere in an infinite volume of fluid accelerates alternately in the
# positive and negative x-direction, generating positive pressure ahead of it,
# negative pressure behind. Time-dependent simulation. Described in [1].

# ## References

# [1] Krysl P, Hawkins AD, Schilt C, Cranford TW (2012): Angular Oscillation of
# Solid Scatterers in Response to Progressive Planar Acoustic Waves: Do Fish
# Otoliths Rock?. PLOS ONE 7(8): e42591. https://doi.org/10.1371/journal.pone.0042591

# ![](sphere_dipole.png)

# ## Goals

# - Show how to generate hexahedral mesh, mirroring and merging together parts.
# - Execute transient simulation by the trapezoidal-rule time stepping of [1].

##
# ## Definitions

# The finite element code relies on the basic functionality implemented in this
# package.
using FinEtools
# The acoustics functionality is brought in:
using FinEtoolsAcoustics

# We shall need some facilities from the linear algebra package
using LinearAlgebra

# The properties correspond roughly to air.
rho = 1.21*phun("kg/m^3");# mass density
c  = 343.0*phun("m/s");# sound speed
bulk =  c^2*rho;
a_amplitude=1.0*phun("mm/s^2");# amplitude of the  acceleration of the sphere
R = 50.0*phun("mm"); # radius of the interior sphere
Ro = 8*R # radius of the external sphere, where the fluid doing this truncated
P_amplitude = R*rho*a_amplitude; # pressure amplitude
frequency = 1200.; # # frequency of the incident wave, Hz
omega = 2*pi*frequency;
dt = 1.0/frequency/20;
tfinal = 90*dt;
nsteps = round(tfinal/dt)+1;

# This dictates the parameters of the mesh.
nref = 2;
nlayers = 40;
tolerance  =  Ro/(nlayers)/100

# Hexahedral mesh: Generate 1/8 of the scattering sphere. This sphere is not
# actually part of the model, it just serves as an intermediate step from which
# the volume of the fluid will be produced.
fens, fes  =  H8sphere(R, nref);
# No mirror the 1/8 of the sphere to build one fourth of the complete
# three-dimensional geometry. In other words we are using two symmetric planes
# here. Note that in order to mirror the hexahedral elements correctly, we need
# to specify the renumbering rule for the connectivity.
fens1, fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0], renumb =  r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]);
# Now merge the two pieces (1/8) together.
fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
fes = cat(newfes1,fes2)

# Derive the finite element set for the boundary of the target.
bfes  =  meshboundary(fes)

# To identify the outer spherical boundary we'll use this function to calculate
# its normal.
function dout(xyz)
    return xyz/norm(xyz)
end
# Select the part of the boundary of the target facing towards infinity.
louter = selectelem(fens, bfes, facing = true, direction = dout)
outer_fes = subset(bfes, louter);

# The surface of the spherical target is now extruded into multiple layers.
fens,fes  = H8extrudeQ4(fens, outer_fes, nlayers, (xyz, layer)->(R+layer/nlayers*(Ro-R))*xyz/norm(xyz));
# In case there are any unconnected nodes, remove them, and renumber the elements.
connected = findunconnnodes(fens, fes);
fens, new_numbering = compactnodes(fens, connected);
fess = renumberconn!(fes, new_numbering);

# The geometry and the solution (pressure) fields.
geom  =  NodalField(fens.xyz)
P  =  NodalField(zeros(FCplxFlt,size(fens.xyz,1),1))

# We need to identify the surface of the scattering sphere, and also the outside
# sphere where the fluid immediately surrounding the sphere is adjacent to the
# infinite extent of the fluid.
bfes  =  meshboundary(fes)

# Find the two spherical surfaces, one can be located by the distance criterion
# (from the center), and the other one by facing to infinity (away from the
# center).
linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
louter = selectelem(fens, bfes, facing = true, direction = dout)

##
# ## Visualize the geometry

# Export three VTK files: one for the interior of the fluid, and two for the boundaries.
# The boundaries on the symmetry planes are taken into account implicitly.
File  =   "Sphere-fluid-interior.vtk"
vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.VTK.H8)
File  =   "Sphere-fluid-inner-boundary.vtk"
vtkexportmesh(File, subset(bfes, linner).conn, geom.values, FinEtools.MeshExportModule.VTK.Q4)
File  =   "Sphere-fluid-outer-boundary.vtk"
vtkexportmesh(File, subset(bfes, louter).conn, geom.values, FinEtools.MeshExportModule.VTK.Q4)

# If we have paraview, we can run it now.
@async run(`"paraview.exe" $File`)

##
# ## Set up the discrete model

# Number the degrees of freedom in the pressure field.
numberdofs!(P)

# Create the finite element machine for the fluid.
material = MatAcoustFluid(bulk,rho)
femm  =  FEMMAcoust(IntegDomain(fes, GaussRule(3, 2)), material)
# Use the machine calculate the acoustic stiffness and mass matrices.
S  =  acousticstiffness(femm, geom, P);
C  =  acousticmass(femm, geom, P);

# Set up finite element machine needed for the absorbing boundary conditions.
abcfemm  =  FEMMAcoustSurf(IntegDomain(subset(bfes, louter), GaussRule(2, 2)), material)
# This is the "damping" mmatrix for the ABC.
D  =  acousticABC(abcfemm, geom, P);

# The sphere at the center is rigid and moves in a prescribed, harmonic,
# fashion. That generates loading onto the fluid. We expect the loading to be
# in the form of a dipole.
function dipole(dpdn, xyz, J, label, t)
    n = cross(J[:,1],J[:,2]);
    n = vec(n/norm(n));
    dpdn[1] = -rho*a_amplitude*sin(omega*t)*n[1]
end
# In order to evaluate the distributed pressure flux loading we need a surface
# acoustic finite element model machine. We shall set it up for the interior
# spherical surface.
dipfemm  =  FEMMAcoustSurf(IntegDomain(subset(bfes, linner), GaussRule(2, 2)), material)

##
# ## Time stepping

# Solve the transient acoustics equations. Refer to [1] for details of the
# formulation.

# The loop executes inside this local scope.
P1 = let
    P0 = deepcopy(P)
    P0.values[:] .= 0.0; # initially all pressure is zero
    vP0 = gathersysvec(P0);
    vQ0 = zeros(eltype(vP0), size(vP0));
    # The `P1` field will be the output of this computation: the final value of
    # the pressure field
    P1 = deepcopy(P0);

    t = 0.0; # Initial time
    # Compute the initial load due to the pressure gradient on the surface of
    # the moving sphere.
    fi  =  ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
    La0 = distribloads(dipfemm, geom, P1, fi, 2);
    # This is the coefficient matrix that needs to be used in the solves. We are
    # not being very careful here to save on computation: it might be best to
    # factorize this matrix, and then use backward and forward solves inside
    # the loop.
    A = (2.0/dt)*S + D + (dt/2.)*C;

    step = 0;
    while t <= tfinal
        step = step  + 1;
        println("Time $t ($(step)/$(round(tfinal/dt)+1))")
        t = t+dt;
        # Compute the current load due to the pressure gradient on the surface of
        # the moving sphere.
        fi  = ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
        La1 = distribloads(dipfemm, geom, P1, fi, 2);
        # Solve for the rate of the pressure
        vQ1 = A\((2/dt)*(S*vQ0)-D*vQ0-C*(2*vP0+(dt/2)*vQ0)+La0+La1);
        # Update the value of the pressure
        vP1 = vP0 + (dt/2)*(vQ0+vQ1);
        # Swap variables for next step
        vP0 = vP1;
        vQ0 = vQ1;
        P1 = scattersysvec!(P1, vec(vP1));
        P0 = P1;
        La0 = La1;
    end

    P1 # Return the final pressure
end

##
# ## Visualization

File  =   "sphere_dipole_P1.vtk"
vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.VTK.H8; scalars = [( "realP", real.(P1.values)),])
@async run(`"paraview.exe" $File`)

true


text:
using Gadfly, Base.Dates

set_default_plot_size(6inch, 3inch)

dates = Array{Date}(40)
dates[1] = today()
for i=2:length(dates)
    dates[i] = today()+Day(i)
end

plot(x=dates,y=sort(rand(40)*40),Geom.bar)


text:
function nlp_mi_004_012(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)

    # Test Goals:
    # - feasible model
    # - no objective
    # - discrete and non-discrete variables

    model = Model(optimizer)

    @variable(model, -1 <= x <= 1)
    @variable(model, y)
    @variable(model, z, Int)

    @NLconstraint(model, x^2 + y^2 + z^2 <= 10)
    @NLconstraint(model, -1.2 * x - y <= z / 1.35)
    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    # 0.0 as it doesn't have an objective
    return check_objective(model, 0.0, tol = objective_tol)
end