text:
# We thank QuTiP developers for sharing the example on Quantum MCWF method,
# this is a direct translation of QuTiP version to solving in Julia using
# QuDynamics.

using QuBase
using QuDynamics

# Parameters
wc = 1.0 * 2 * pi   # cavity frequency
wa = 1.0 * 2 * pi   # atom frequency
g  = 0.25 * 2 * pi  # coupling strength
kappa = 0.05        # cavity dissipation rate
gamma = 0.15        # atom dissipation rate
use_rwa = true
tlist = linspace(0., 10., 200)

# System construction

# Hamiltonian
idc = QuArray(eye(10))
ida = QuArray(eye(2))
a  = tensor(ida, lowerop(10))
sm = tensor(lowerop(2), idc)
if use_rwa
    # use the rotating wave approxiation
    hamiltonian = wc * a' * a + wa * sm' * sm + g * (a' * sm + a * sm')
else
    hamiltonian = wc * a' * a + wa * sm' * sm + g * (a' + a) * (sm + sm')
end

# initial state
init_state = complex(tensor(statevec(1, FiniteBasis(2)), statevec(6, FiniteBasis(10))))
init_state_dm = complex(init_state*init_state')

# collapse operators
c_ops = [sqrt(0.1) * a]

# rho being constructed to store density matrix at every time step
rho = Array(typeof(init_state_dm), length(tlist)-1)
qumcwfen = QuMCWFEnsemble(complex(init_state), 1000)
for i=1:length(tlist)-1
    rho[i] = complex(zeros(init_state_dm))
end

# Using MCWF and extracting the density matrices
for psi0 in qumcwfen
    i = 1
    mcwf = QuMCWF()
    for (t,psi) in QuPropagator(hamiltonian, c_ops, psi0, tlist, mcwf)
        rho[i] = rho[i] + (psi*psi')/length(qumcwfen)/norm(psi)^2
        i = i + 1
    end
end

# Calculating the expectation of the observable a'*a and sm'*sm at every time step
for i in 1:length(tlist)-1
    println("t = $(tlist[i+1]) : <a'*a> = $(trace(rho[i]*a'*a)), <sm'*sm> = $(trace(rho[i]*sm'* sm))")
end


text:
__precompile__(true)

"""
A general purpose Gaussian Process implementation for Julia.
"""
module GPs

include("general.jl")

#--- Other submodules

include("Metrics.jl")
include("Kernels.jl")
include("DataSets.jl")
include("GaussianLikelihoods.jl")

using .Metrics
using .Kernels
using .DataSets
using .GaussianLikelihoods

export GP, TrainingSet, SqrExponential, Euclidean

end # module

ext: jl


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

JLLWrappers.@generate_main_file_header("primecount")
JLLWrappers.@generate_main_file("primecount", UUID("ba3b429d-c409-5a27-9e45-21aab9233266"))
end  # module primecount_jll


text:
module TestTriangle

using Test
using Triangle

include("test_Triangle_NativeInterface_trioptions.jl")
include("test_Triangle.jl")

end


text:
@testset "invalidate native code cache" begin
    # invalidate native code cache in a system image if it has not been analyzed by JET
    # yes this slows down anlaysis for sure, but otherwise JET will miss obvious errors like below
    let
        interp, frame = analyze_call((Nothing,)) do a
            a.field
        end
        @test length(interp.reports) === 1
    end

    # invalidation from deeper call site can refresh JET analysis
    let
        # NOTE: branching on https://github.com/JuliaLang/julia/pull/38830
        symarg = last(first(methods(Base.show_sym)).sig.parameters) === Symbol ?
                 :(sym::Symbol) :
                 :(sym)

        l1, l2, l3 = @freshexec begin
            # ensure we start with this "errorneous" `show_sym`
            @eval Base begin
                function show_sym(io::IO, $(symarg); allow_macroname=false)
                    if is_valid_identifier(sym)
                        print(io, sym)
                    elseif allow_macroname && (sym_str = string(sym); startswith(sym_str, '@'))
                        print(io, '@')
                        show_sym(io, sym_str[2:end]) # NOTE: `sym_str[2:end]` here is errorneous
                    else
                        print(io, "var", repr(string(sym)))
                    end
                end
            end

            # should have error reported
            interp1, = @analyze_call println(QuoteNode(nothing))

            # should invoke invalidation in the deeper call site of `println(::QuoteNode)`
            @eval Base begin
                function show_sym(io::IO, $(symarg); allow_macroname=false)
                    if is_valid_identifier(sym)
                        print(io, sym)
                    elseif allow_macroname && (sym_str = string(sym); startswith(sym_str, '@'))
                        print(io, '@')
                        show_sym(io, Symbol(sym_str[2:end]))
                    else
                        print(io, "var", repr(string(sym)))
                    end
                end
            end

            # now we shouldn't have reports
            interp2, = @analyze_call println(QuoteNode(nothing))

            # again, invoke invalidation
            @eval Base begin
                function show_sym(io::IO, $(symarg); allow_macroname=false)
                    if is_valid_identifier(sym)
                        print(io, sym)
                    elseif allow_macroname && (sym_str = string(sym); startswith(sym_str, '@'))
                        print(io, '@')
                        show_sym(io, sym_str[2:end])
                    else
                        print(io, "var", repr(string(sym)))
                    end
                end
            end

            # now we should have reports, again
            interp3, = @analyze_call println(QuoteNode(nothing))

            length(interp1.reports), length(interp2.reports), length(interp3.reports) # return
        end

        @test l1 > 0
        @test l2 == 0
        @test l3 == l1
    end
end

@testset "integrate with global code cache" begin
    # analysis for `sum(::String)` is already cached, `sum′` and `sum′′` should use it
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            sum′(s) = sum(s)
            sum′′(s) = sum′(s)
            $analyze_call() do
                sum′′("julia")
            end
        end)
        test_sum_over_string(interp)
    end

    # incremental setup
    let
        m = gen_virtual_module()

        interp, frame = Core.eval(m, quote
            $analyze_call() do
                sum("julia")
            end
        end)
        test_sum_over_string(interp)

        interp, frame = Core.eval(m, quote
            sum′(s) = sum(s)
            $analyze_call() do
                sum′("julia")
            end
        end)
        test_sum_over_string(interp)

        interp, frame = Core.eval(m, quote
            sum′′(s) = sum′(s)
            $analyze_call() do
                sum′′("julia")
            end
        end)
        test_sum_over_string(interp)
    end

    # should not error for virtual stacktrace traversing with a frame for inner constructor
    # https://github.com/aviatesk/JET.jl/pull/69
    let
        res = @analyze_toplevel begin
            struct Foo end
            println(Foo())
        end
        @test isempty(res.inference_error_reports)
    end
end

@testset "integrate with local code cache" begin
    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            struct Foo{T}
                bar::T
            end
            $analyze_call((Foo{Int},)) do foo
                foo.baz # typo
            end
        end)

        @test !isempty(interp.reports)
        @test !isempty(interp.cache)
        @test any(interp.cache) do analysis_result
            analysis_result.argtypes==Any[CC.Const(getproperty),m.Foo{Int},CC.Const(:baz)]
        end
    end

    let
        m = gen_virtual_module()
        interp, frame = Core.eval(m, quote
            struct Foo{T}
                bar::T
            end
            getter(foo, prop) = getproperty(foo, prop)
            $analyze_call((Foo{Int}, Bool)) do foo, cond
                getter(foo, :bar)
                cond ? getter(foo, :baz) : getter(foo, :qux) # non-deterministic typos
            end
        end)

        # there should be local cache for each errorneous constant analysis
        @test !isempty(interp.reports)
        @test !isempty(interp.cache)
        @test any(interp.cache) do analysis_result
            analysis_result.argtypes==Any[CC.Const(m.getter),m.Foo{Int},CC.Const(:baz)]
        end
        @test any(interp.cache) do analysis_result
            analysis_result.argtypes==Any[CC.Const(m.getter),m.Foo{Int},CC.Const(:qux)]
        end
    end
end


text:
abstract type AbstractMicrophysics <: AbstractPhysicsComponent end

Base.@kwdef struct ZeroMomentMicrophysics <: AbstractMicrophysics end

@inline calc_component!(source, ::Nothing, state, _...) = nothing

@inline function calc_component!(source, ::ZeroMomentMicrophysics, state, aux, physics)
    ρ = state.ρ
    ρq = state.ρq
    Φ = aux.Φ
    eos = physics.eos
    τ = physics.parameters.τ_precip
    T_0 = physics.parameters.T_0
    cv_l = physics.parameters.cv_l

    # we need the saturation excess in order to calculate the
    # source terms for prognostic variables
    T = calc_air_temperature(eos, state, aux, physics.parameters)
    qᵥₛ = calc_saturation_specific_humidity(ρ, T, physics.parameters)
    ρS = max(0, ρq - ρ * qᵥₛ) # saturation excess
    Iₗ = cv_l * (T - T_0) # liquid internal energy

    # source terms are proportional to the saturation excess
    source.ρ -= ρS / τ
    source.ρe -= (Iₗ + Φ) * ρS / τ
    source.ρq -= ρS / τ
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
using AdFem
m = 5
n = 5
h = 1/n

# F(u + Delta u, v) = F_u Delta u + F

# F(u , v + Delta v) = F + F_v Delta v

function eval_f(S)
    u, v = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    # ux = ones(4*m*n)

    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    return f1 + f2
end

function eval_df(S)
    u, v = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = compute_fem_mass_matrix1(ux, m, n, h)
    M2 = run(sess, compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    Fv = compute_fem_mass_matrix1(uy, m, n, h)
    return M1+M2, Fv
end

sess = Session()
S = rand(2(m+1)*(n+1))
function fJ(x)
    global S
    S[(m+1)*(n+1)+1:end] = x
    f = eval_f(S)
    _, df = eval_df(S)
    return f, df
end

function fJ2(x)
    global S
    S[1:(m+1)*(n+1)] = x
    f = eval_f(S)
    df, _ = eval_df(S)
    return f, df
end

function ADCME.:test_jacobian(f::Function, x0::Array{Float64}; scale::Float64 = 1.0)
    v0 = rand(Float64,size(x0))
    γs = scale ./10 .^(1:5)
    err2 = []
    err1 = []
    f0, J = f(x0)
    for i = 1:5
        f1, _ = f(x0+γs[i]*v0)
        push!(err1, norm(f1-f0))
        @show f1, f0, 2γs[i]*J*v0
        push!(err2, norm(f1-f0-γs[i]*J*v0))
        # push!(err2, norm((f1-f2)/(2γs[i])-J*v0))
        # #@show "test ", f1, f2, f1-f2
    end
    close("all")
    loglog(γs, err2, label="Automatic Differentiation")
    loglog(γs, err1, label="Finite Difference")
    loglog(γs, γs.^2 * 0.5*abs(err2[1])/γs[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(γs, γs * 0.5*abs(err1[1])/γs[1], "--",label="\$\\mathcal{O}(\\gamma)\$")
    plt.gca().invert_xaxis()
    legend()
    println("Finite difference: $err1")
    println("Automatic differentiation: $err2")
    return err1, err2
end

u0 = rand((m+1)*(n+1))
fJ(u0)

close("all")
test_jacobian(fJ2, u0)


text:
import LinearAlgebra
using LinearAlgebra: norm, dot, I, UniformScaling, Symmetric, eigvals, eigen
using Printf

export dgmres, dlanczos, darnoldi

"""
`dirder(x,w,f,f0)`

Finite difference directional derivative, approximate f'(x) w

### Inputs
* x, w = point and direction
* f = function
* f0 = f(x), in nonlinear iterations f(x) has usually been computed
                before the call to dirder

### Output
* `z`: approximation to `f'(x) w`

based on a code by C. T. Kelley, April 1, 2003
"""
function dirder(x, w, f, f0, h=1e-7)
   # initialise difference increment
   epsnew = h
   n = length(x)
   # scale the step
   if norm(w) == 0
      return zeros(n)
   end
   # Now scale the difference increment.
   xs = dot(x, w) / norm(w)
   if xs != 0.0
      epsnew *= max(abs(xs), 1.0) * sign(xs)
   end
   epsnew /= norm(w)
   f1 = f(x + epsnew * w)
   return (f1 - f0) / epsnew
end

"""
`dirderinf(f, f0, xc, v, h)`:
directional derivative operator to mimic the ∂f(xc) * v action
via (f(x+h'*v) - f(x)) / h' where h' is rescaled against ||v||_∞
"""
function dirderinf(f, f0, xc, v, h=1e-7)
   h = h / norm(v, Inf)
   return (f(xc + h * v) - f0) / h
end

"""
`function givapp(c, s, vin, k) -> vrot`

Apply a sequence of k Givens rotations, used within gmres codes.
"""
function givapp(c, s, vin, k)
   vrot = copy(vin)                   # TODO: can we do this in-place?
   for i = 1:k
      w1 = c[i] * vrot[i] - s[i] * vrot[i+1]
      # Here's a modest change that makes the code work in complex
      # arithmetic. Thanks to Howard Elman for this.
      # w2 = s(i)*vrot(i)+c(i)*vrot(i+1);
      w2 = s[i] * vrot[i] + c[i]' * vrot[i+1]
      vrot[i:i+1] = [w1, w2]
   end
   return vrot
end

"""
`function dgmres(f0, f, xc, errtol, kmax, [reorth, xinit]) -> x, error, total_iters`

GMRES linear equation solver for use in Newton-GMRES solver

### Input Parameters

* f0 : function at current point, f(xc)
* f : nonlinear function the format for f is  function fx = f(x)
              Note that for Newton-GMRES we incorporate any
              preconditioning into the function routine
* xc : current point
* errtol : relative residual reduction factor
* kmax : max number of iterations
* reorth : reorthogonalization method
      - 1 -- Brown/Hindmarsh condition (default)
      - 2 -- Never reorthogonalize (not recommended)
      - 3 -- Always reorthogonalize (not cheap!)
* xinit : initial iterate. xinit = zeros(f0) is the default. This
              is a reasonable choice unless restarted GMRES
              will be used as the linear solver.

### Output Parameters

* x : solution
* error : vector of residual norms for the history of the iteration
* total_iters : number of iterations
"""
function dgmres(f0, f, xc, errtol, kmax, reorth = 1, x = zeros(length(f0)); hfd = 1e-7)
   # The right side of the linear equation for the step is -f0.
   b = - f0
   n = length(b)
   r = - dirder(xc, x, f, f0, hfd) - f0

   h = zeros(kmax+1, kmax+1)
   v = zeros(n, kmax+1)
   c = zeros(kmax+1)
   s = zeros(kmax+1)
   rho = norm(r)
   g = rho * Matrix(1.0*I, kmax+1, 1)
   errtol = errtol * norm(b)
   error = Float64[]

   # Test for termination on entry.
   push!(error, rho)
   total_iters = 0
   if (rho < errtol)
      # early termination
      return x, error, total_iters
   end

   v[:,1] = r / rho
   beta = rho
   k = 0

   # GMRES iteration
   while rho > errtol && k < kmax
      k = k+1;

      # Call directional derivative function.
      v[:, k+1] = dirder(xc, v[:,k], f, f0, hfd)
      normav = norm(v[:,k+1])

      # Modified Gram-Schmidt
      for j = 1:k
         h[j, k] = dot(v[:,j], v[:,k+1])
         v[:, k+1] = v[:,k+1] - h[j,k] * v[:,j]
      end
      h[k+1,k] = norm(v[:,k+1])
      normav2 = h[k+1,k]

      # Reorthogonalize?
      if  ((reorth == 1) && (normav + .001 * normav2 == normav)) || (reorth ==  3)
         for j = 1:k
            hr = dot(v[:,j], v[:,k+1])
            h[j,k] = h[j,k] + hr
            v[:,k+1] = v[:,k+1] - hr * v[:,j]
         end
         h[k+1, k] = norm(v[:,k+1])
      end

      # Watch out for happy breakdown.
      if h[k+1,k] != 0
         v[:,k+1] = v[:,k+1] / h[k+1,k]
      end

      # Form and store the information for the new Givens rotation.
      if k > 1
         h[1:k,k] = givapp(c[1:k-1], s[1:k-1], h[1:k,k], k-1)
      end

      # Don't divide by zero if solution has  been found.
      nu = norm(h[k:k+1,k])
      if nu != 0
         # c[k] = h[k,k]/nu
         c[k] = conj(h[k,k] / nu)
         s[k] = -h[k+1,k] / nu;
         h[k,k] = c[k] * h[k,k] - s[k] * h[k+1,k]
         h[k+1,k] = 0
         g[k:k+1] = givapp(c[k], s[k], g[k:k+1], 1)
      end

      # Update the residual norm.
      rho = abs(g[k+1])
      push!(error, rho)
   end # end of the main while loop

   # At this point either k > kmax or rho < errtol.
   # It's time to compute x and leave.
   y = h[1:k,1:k] \ g[1:k]
   total_iters = k
   x = x + v[1:n, 1:k] * y

   return x, error, total_iters, total_iters
end

"""
`pushcol(V::Matrix, v::Vector) -> V`:
append the vector `v` as a column to `V`
"""
function pushcol(V::Matrix, v::Vector)
   rows, cols = size(V)
   V = V[:]
   append!(V, v)
   return reshape(V, rows, cols+1)
end

"""
`sorted_eig(A::Union{SymTridiagonal, Symmetric})`:
diagonalise A, and return sorted eigenvalues and analogously sorted eigenvectors
"""
function sorted_eig(A)
   D, Q = eigen(A)
   I = sortperm(real(D))
   return D[I], Q[:, I]
end

LinearAlgebra.dot(x, A::UniformScaling, y) = A.λ * dot(x,y)
LinearAlgebra.dot(x, A::AbstractMatrix, y) = dot(x, A*y)
LinearAlgebra.norm(P, x) = sqrt(abs(dot(x, P*x)))
dualnorm(P, f) = sqrt(abs(dot(f, P \ f)))
# nkdualnorm(P, f) = dualnorm(P, f)  # TODO: revisit this!!!
nkdualnorm(P, f) = norm(f)
# There seems to be no consistent benefit to either!

"""
`orthogonalise(w, V::Matrix, P=I) -> w'`:
Assuming V is P-orthogonal, P-orthogonalise w to V and return; if
w is already in the span of V, then a random vector is orthogonalised
and returned instead. The returned vector is normalised.
"""
function orthogonalise(w, V::Matrix, P=I; maxiter = 10, ORTHTOL=1e-10)
   @assert size(V,2) < size(V,1)
   for n = 1:maxiter
      w /= norm(P, w)
      # make w orthogonal to all columns of V
      for i = 1:size(V,2)
         w -= dot(w, P, V[:,i]) * V[:,i]
      end
      # test that the resulting w is truly orthogonal:
      nrmw = norm(P, w)
      if nrmw < 1e-15 && n == maxiter; break; end
      if nrmw < 1e-15
         w = P \ (rand(length(w)) - 0.5)
         w /= norm(P, w)
      else
         w /= norm(P, w)
         if norm( (P*w)' * V, Inf ) < ORTHTOL
            return w
         end
      end
   end
   warn("""`orthogonalise` did not terminate succesfully, most likely
      the returned vector is not orthogonal to V""")
   return w
end

"""
`appendkrylov(V, AxV, Y, v, Hmul, P)`:
appends [V v], [AxV A*v], [Y P \\ A*v]
"""
function appendkrylov(V, AxV, Y, v, Hmul, P)
   V   = pushcol(V,   v)
   AxV = pushcol(AxV, Hmul(v))
   Y   = pushcol(Y,   P \ AxV[:, end])
   return V, AxV, Y
end

"""
`immutable LanczosMatrix`: (todo: write documentation)
"""
struct LanczosMatrix
   V::Matrix{Float64}
   D::Vector{Float64}
   P
end

Base.length(K::LanczosMatrix) = length(K.D)
Base.getindex(K::LanczosMatrix, i) = K.D[i], K.V[:,i]
import Base: *, \
import LinearAlgebra: mul!
*(K::LanczosMatrix, v) = K.P * (K.V * (K.D .* (K.V' * (K.P * v))))
mul!(out::Array{Float64,1}, L::LanczosMatrix, x::Array{Float64,1}) = copy!(out, L * x)
\(K::LanczosMatrix, f) = K.V * ( K.D.^(-1) .* (K.V' * f) )
# rank
# isnewton

"spectrum transformation for `dlanczos` for minimisation"
minim(D) = abs.(D)
"spectrum transformation for `dlanczos` for index-1 saddles"
index1(D) = [-abs(D[1]); abs.(D[2:end])]
"spectrum transformation for `dlanczos` for index-2 saddles"
index2(D) = [-abs.(D[1:2]); abs.(D[3:end])]

"spectrum transformation for general order saddle"
indexp(D, p::Integer) = [-abs.(D[1:p]); abs.(D[p+1:end])]

# TODO:
#   * look into re-orthogonalising
#   * look into getting rid of AxV
#   * this code does not make proper use of the lanczos recursion
#     at all; either switch to an arnoldi method, or clean it up
#   * `isnewton` could be incorporated into LanczosMatrix

"""
dlanczos(f0, f, xc, errtol, kmax, transform; kwargs...)
      -> x, G, numf, success, isnewton

Given E : ℝᴺ → ℝ with hessian H = ∇²E(x) and b ∈ ℝᴺ  `dlanczos` computes an
approximate solution `u` to the system
   g(H) u = b
where g(H) ∈ ℝᴺˣᴺ is a transformation of J = Q D Q'  of the form
   g(H) = Q g(D) Q'
Here, D = diag(λ₁, λ₂, …) is the diagonal matrix of ordered eigenvalues
(λ₁ being the smallest). More on the transformation g see below.

`dlanczos` first uses a (possibly preconditioned) (block-)lanczos method to compute
a reasonable approximation to H in the form H = P V T V' P. Then it diagonalises
T = Q D Q', replaces D with g(D) as above. This gives an approximation G to g(H).
The method terminates if `pinv(G) b` yields a solution of sufficiently high accuracy
(residual).

In the lanzcos iterations the matrix vector product H * u
is replaced with the finite-difference operation (∇E(xc + h u) - ∇E(xc))/h;
see also `dirder` and `dirderinf`.

## Required Parameters

* `f0` : ∇E(xc)
* `f` : function to evaluate ∇E
* `xc` : current point
* `b` : right-hand side in the equation
* `errtol` : max-norm residual tolerance for termination
* `kmax` : maximum number of lanzcos iterations
* `transform` : the g transformation, default: id

## KW parameters

* `P` : preconditioner (default: I;  minimally needs to define `*` and `\\` )
* `V0` : initial subspace (default: P \\ b)
* `eigatol`, `eigrtol`: tolerance on the first eigenvalue
* `debug`: show debug information (true/false)
* `hfd` : finite-difference parameter
* `dirder` : function to compute the directional derivative; see
         `dirderinf` for format
* `ORTHTOL` : orthogonality tolerance parameter
* `Hmul` : the finite-difference operation can in principle be replaced by
         and arbitrary operator specifying the matrix-vector product

## Returns

* `x` : approximate solution
* `G` : of type `LanczosMatrix`, to extract information about the computed operator
* `numf` : number of f evaluations
* `success` : true if termination criteria are satisfied, false if kmax is reached
* `isnewton` : true if the step if the hessian spectrum was left unmodified

## Further Notes

* Acknowledgement: This code is inspired by C T Kelley's Matlab Newton-Krylov solver `nsoli.m`;
see also `nsoli` and `dgmres` which are directly ported with permission from `nsoli.m`.
* Termination criterion: If g ≠ id then G * u - b is not in fact the residual
of the equation we are trying to solve but only an approximation to that residual.
Therefore the termination criterion is only approximately satisfied.
"""
function dlanczos(f0, f, xc, b, errtol, kmax, transform = identity;
                  P = I, V0 = P \ b,
                  eigatol = 1e-1, eigrtol = 1e-1,
                  debug = false, hfd = 1e-7,
                  dirder = dirderinf,
                  Hmul = z -> dirder(f, f0, xc, z, hfd),
                  ORTHTOL = 1e-12, nevals=1 )

   # make sure V0 is given in the correct format (Matrix)
   if isa(V0, Vector)
      V0 = reshape(V0, length(V0), 1)
   end

   # initialise some variables
   p = size(V0, 2)     # block size
   d = length(f0)      # problem dimension
   @assert kmax <= d
   numf = 0            # count f evaluations
   isnewton = false    # remember whether the output is a newton direction
   success = false     # flag for termination

   # initialise Krylov subspace and more
   V = zeros(d,0)      # store the Krylov basis
   AxV = zeros(d, 0)   # store A vⱼ
   Y = zeros(d, 0)     # store P \ A vⱼ
   Q = Matrix{Float64}(undef, 0,0)
   E = Vector{Float64}(undef, 0)

   # initialise Krylov subspace; TODO: detect if a vector is in the span of previous ones
   for j = 1:size(V0, 2)
      vj = orthogonalise(V0[:,j], V, P)
      V, AxV, Y = appendkrylov(V, AxV, Y, vj, Hmul, P)
      numf += 1
   end

   # prepare for the Block-Lanczos loop
   j = 1     # the index of the krylov vector which will next be multiplied by J = A
   x = zeros(d)
   λ = sort(eigvals(Symmetric(V'*AxV)))[1:nevals]
   λ_old = Vector{Float64}[]

   if debug
      @printf("      numf     λ    err_λ     |Ax-b|/|b| \n")
   end

   # start the block-lanczos loop; when we have kmax v-vectors we stop
   while size(V, 2) <= kmax

      # At this point we have V, AxV and Y = P \ AxV  available, so we can
      # now solve the projected linear system and eigenvalue problem
      #
      # there is probably an elegant way to assemble the projected
      # linear system, but for now we just do it brute-force:
      #    A vj = P V T V' P vj
      #    yj = P \ A qj = V T  V' P vj
      #    V' P Y = V' P V T V' P V = T
      n = size(V, 2)
      T = Symmetric(V' * AxV)       # T = V' * P * Y = V' * AV
      #  TODO: could replace above with with [dot(V[:,a], P, Y[:,b] for a=1:n, b=1:n]
      #        to avoid storing AxV altogether

      # make the specified transformation
      D, Q = sorted_eig(T)
      E = transform(D)
      # residual estimate for the old x    >>> TODO: should we switch to P^{-1}-norm?
      res = nkdualnorm(P,  P * (V * (Q * (E .* (Q' * (V' * (P * x)))))) - b )
      # new x and λ (remember the old)
      g = Q * (E .\ (Q' * (V' * b)))
      x, x_old = V * g, x
      push!(λ_old, λ)
      λ = D[1:nevals]
      # if E == D then g(H) = H hence we can estimate the *actual* and *current* residual
      if (isnewton = (norm(E - D, Inf) < 1e-7) )
         res = nkdualnorm(P, AxV * g - b )    # TODO: should we switch to P^{-1}-norm?
      end
      # check for termination
      err_λ = maximum(norm(λ - λ_old[i], Inf)
                      for i = max(1,length(λ)-p+1):length(λ))
      if debug
         @printf("      %d   %.2f   %.2e    %.2e \n", numf, λ[1], err_λ, res/norm(b))
      end

      # CHECK FOR SUCCESFUL TERMINATION
      #  * if size(V,2) == d then we have assembled the entire matrix
      if (size(V,2) == d) ||  (
         (res < errtol) && ((λ == E[1:nevals]) || (err_λ < eigatol + eigrtol * abs(λ[1])))  )
         success = true
         break
      end

      # CHECK FOR UNSUCCESFUL TERMINATION
      if size(V,2) == kmax
         success = false
         break
      end

      # add the next Krylov vector
      if j <= size(V,2)   # we have Krylov vectors left to cycle through (99.9% of the time)
         w = orthogonalise(Y[:, j], V, P)
         V, AxV, Y = appendkrylov(V, AxV, Y, w, Hmul, P)
         numf += 1
         j += 1
      else
         error("how did we get here?")
      end
      # else    # no vectors left, try to add a random vector (how can this happen???)
      #    warn("no vectors left; had to add a random vector")
      #    w = orthogonalise(P \ rand(d), V, P)
      #    V, AxV, Y = appendkrylov(V, AxV, Y, w, Hmul, P)
      #    numf += 1
      # end
   end
   # if we are here it means that kmax is reached, i.e. we terminate with
   # warning or error >>> TODO: return to how to best handle this?!
   # warn("`dcg_index1` did not converge within kmax = $(kmax) iterations")
   return x, LanczosMatrix(V * Q, E, P), numf, success, isnewton
end

# """
# `immutable KrylovMatrix`: (todo: write documentation)
#
# We represent a matrix projected to a Krylov subspace as
#
# P * V * X * diagm(D) * X⁻¹ * V' * P
#
# and its inverse as
#
# V * X * diagm(D.⁻¹) * X⁻¹ * V'
# """
# immutable KrylovMatrix{T}
#    V::Matrix{T}   # orthogonal component
#    X::Matrix{T}   # non-orthogonal component
#    D::Vector{T}   # diagonal component
#    P              # preconditioner (we don't know anything about this object)
# end
#
# Base.length(K::KrylovMatrix) = length(K.D)
# Base.rank(K::KrylovMatrix) = length(K)
# Base.getindex(K::KrylovMatrix, i) = K.D[i], K.V[:,i]
# import Base: *, \, A_mul_B!
# *(K::LanczosMatrix, v) = K.P * (K.V * (K. X * (K.D .* (K.X \ (K.V' * (K.P * v))))))
# A_mul_B!(out::Array{Float64,1}, L::LanczosMatrix, x::Array{Float64,1}) = copy!(out, L * x)
# \(K::LanczosMatrix, f) = K.V * ( K.D.^(-1) .* (K.V' * f) )

"""
`stabilise(D::Vector) = D .* sign.(real.(D))`

a Lipschitz-continuous spectrum transformation that stabilises the origin
"""
stabilise(D::Vector) = D .* sign.(real.(D))

"""
darnoldi(f0, f, xc, errtol, kmax, transform; kwargs...)
      -> u, G, numf, success, isnewton

Given f : ℝᴺ → ℝᴺ with jacobian J = ∂f(x) and b ∈ ℝᴺ  `darnoldi` computes an
approximate solution `u` to the system
   g(J) u = b
where g(J) ∈ ℝᴺˣᴺ is a transformation of J = X D X⁻¹  of the form
   g(J) = X g(D) X⁻¹
Here, D = diag(λ₁, λ₂, …) is the diagonal matrix of ordered eigenvalues
(λ₁ being the smallest). More on the transformation g see below. Also
see below what happens if J is not diagonalisable.

`darnoldi` first uses a (possibly preconditioned) (block-)arnoldi iteration to compute
an approximation to  J in the form J = P V H V'. Then it diagonalises
H = X D V⁻¹, replaces D with g(D) as above. This gives an approximation G to g(J).
The method terminates if `pinv(G) b` yields a solution of sufficiently high accuracy
(residual).

In the Arnoldi iterations the matrix vector product J * u
is replaced with the finite-difference operation (f(xc + h u) - f(xc))/h;
see also `dirder` and `dirderinf`.

## Required Parameters

* `f0` : f(xc)
* `f` : function to evaluate ∇E
* `xc` : current point
* `b` : right-hand side in the equation
* `errtol` : max-norm residual tolerance for termination
* `kmax` : maximum number of lanzcos iterations
* `transform` : the g transformation, default: id

## KW parameters

* `P` : preconditioner (default: I;  minimally needs to define `*` and `\\` )
* `V0` : initial subspace (default: P \\ b)
* `debug`: show debug information (true/false)
* `hfd` : finite-difference parameter
* `dirder` : function to compute the directional derivative; see
         `dirderinf` for format
* `ORTHTOL` : orthogonality tolerance parameter
* `Hmul` : the finite-difference operation can in principle be replaced by
         an arbitrary operator specifying the matrix-vector product

## Returns

* `x` : approximate solution
/** `G` : of type `LanczosMatrix`, to extract information about the computed operator*/
* `numf` : number of f evaluations
* `success` : true if termination criteria are satisfied, false if kmax is reached
* `isnewton` : true if the hessian spectrum was left unmodified

For Further Notes see `?dlanczos`
"""
function darnoldi( f0, f, xc, b, errtol, kmax, transform = identity;
                   P = I, V0 = P \ b,
                   debug = false, hfd = 1e-7,
                   dirder = dirderinf,
                   Hmul = z -> dirder(f, f0, xc, z, hfd),
                   ORTHTOL = 1e-10 )

   # make sure V0 is given in the correct format (Matrix)
   if isa(V0, Vector)
      V0 = reshape(V0, length(V0), 1)
   end

   # initialise some variables
   p = size(V0, 2)     # block size
   d = length(f0)      # problem dimension
   @assert kmax <= d
   numf = 0            # count f evaluations
   isnewton = false    # remember whether the output is a newton direction
   success = false     # flag for termination

   # initialise Krylov subspace and more
   V = zeros(d,0)      # store the Krylov basis
   AxV = zeros(d, 0)   # store A vⱼ
   Y = zeros(d, 0)     # store P \ A vⱼ
   X = Matrix{ComplexF64}()     # H = X * D * inv(X) becomes
   Dmod = Vector{ComplexF64}()  #    Hmod = X * Dmod * inv(X)

   # initialise Krylov subspace;
   for j = 1:size(V0, 2)
      vj = orthogonalise(V0[:,j], V)
      if dot(vj, V0[:,j]) < ORTHTOL
         warn("the initial subspace matrix V0 does not have full rank")
      end
      V, AxV, Y = appendkrylov(V, AxV, Y, vj, Hmul, P)
      numf += 1
   end

   # prepare for the Arnoldi loop
   j = 1
   x = zeros(d)
   bP = P \ b

   if debug
      @printf("      numf     |Ax-b|/|b| \n")
   end

   # start the block-arnoldi process; when we have kmax v-vectors we stop
   while size(V, 2) <= kmax

      # At this point we have V, AxV and Y = P \ AxV  available, so we can
      # now solve the projected linear system and eigenvalue problem
      n = size(V, 2)
      H = V' * Y
      # make the specified spectrum transformation
      D, X = eigen(H)
      Dmod = transform(D)
      # residual estimate for the old x
      res = norm(  P * (V * (X * (Dmod .* (X \ (V'*x))))) - b )
      # new x (remember the old)
      #     A ≈ P V V' * Y V' = P V H V'     =>     A⁻¹ ≈ V H⁻¹ V' P⁻¹
      xV = X * (Dmod .\ (X \ (V' * bP)))  # coefficients of x in V
      x, x_old = V * xV, x
      # if Dmod == D then this was a newton step
      isnewton = (norm(Dmod - D, Inf) < 1e-7)
      if isnewton
         res = norm(AxV * xV - b)
      end
      if debug
         @printf("      %d    %.2e   %.2e \n", numf, res/norm(b), norm(x-x_old, Inf))
      end

      # CHECK FOR SUCCESFUL TERMINATION
      #  * if size(V,2) == d then we have assembled the entire matrix
      if (size(V,2) == d) || (res < errtol)
         success = true
         break
      end

      # CHECK FOR UNSUCCESFUL TERMINATION
      if size(V,2) == kmax
         success = false
         break
      end

      # add the next Krylov vector
      if j <= size(V,2)   # we have Krylov vectors left to cycle through
         w = orthogonalise(Y[:, j], V)
         V, AxV, Y = appendkrylov(V, AxV, Y, w, Hmul, P)
         numf += 1 # (the call is in appendkrylov)
         j += 1
      else
         error("how did we get here?")
      end
   end
   return x, numf, success, isnewton
end


text:
# Defines the functions used to process the model parameters

# Statistical functions ######################################################
"""
Differential eqn formula for finding densities of each pop at next step
    :param indiv: vector of number of individuals of each group
    (first idx is uninfected hosts; 2:end to each strain)
    :param range: window for each time step
    :param parameters: model parameters
    :returns: differential eqn to solve
"""
function comp_densities(indiv, range, parameters)
    x = indiv[1] # number of uninfected hosts
    y = indiv[2:end] # number of hosts infected with each strain
    regul = (1 .- (sum(indiv)) ./ (range.K))
    dx =
        (range.bx .* x .+ sum(range.ei .* y)) .* regul .- range.ux .* x .-
        range.c .* sum(range.βy .* y) .* x
    dy = range.bi .* y .* regul .- range.ui .* y .+ range.c .* range.βy .* x .* y
    return vcat(dx, dy)
end

"""
Runs the ODE simulation of a population-dynamical model given a set of
parameters. This is done for each strain introduction.
"""
function run_simulation()
    @showprogress for i in 2:length(Y)
        # solve ODE using initial conditions
        prob = ODEProblem(
            comp_densities,
            new_U,
            (windowstart, windowend),
            parameters
        )
        solution = solve(prob; saveat=windowstart:1.0:windowend)
        # set values to 0 if negative
        for t in eachindex(solution.t)
            pop = solution.u[t]
            for i in 1:n_parasites
                if (pop .< 0)[i]
                    pop[i] = 0
                end
            end
            N[:, Int(solution.t[t] + 1)] = pop
        end
        # setup conditions & new parasite for next loop
        global new_U = solution[end]
        new_y = findfirst(x -> x == 0.0, new_U)
        new_U[new_y] = 1.0
        # set limits for next loop
        global windowstart = windowend
        global windowend = windowstart + windowsize
    end
end

"""
Calculates the average mortality across the population
(weighted for noise reduction)
    :returns: vector of mortality values through time
"""
function calculate_mortality()
    return sum(Np[:, 2:end] .* ui'; dims=2) ./ sum(Np[:, 2:end]; dims=2);
end

"""
Calculates the average horizontal transmission (βy) across the population
(weighted for noise reduction)
    :returns: vector of βy values through time
"""
function calculate_horizontal_transmission()
    return sum(Np[:, 2:end] .* βy'; dims=2) ./ sum(Np[:, 2:end]; dims=2);
end

"""
Calculates the average the number of new horizontally-acquired cases from a
single infected host before the primary host dies (H0) across the population
(weighted for noise reduction)
    :returns: vector of H0 values through time
"""
function calculate_H0()
    k = 1
    return c * βy_w_avg ./ ui_w_avg .* k .* (1 - ux / bx);
end

"""
Calculates the vertical basic reproductive ratio (V0) across the population
(weighted for noise reduction)
    :returns: vector of H0 values through time
"""
function calculate_V0()
    bi_avg = sum(Np[:, 2:end] .* bi'; dims=2) ./ sum(Np[:, 2:end]; dims=2);
    return bi_avg .* ux ./ (bx * ui_w_avg)
end

"""
Calculates the virulence across the population (weighted for noise reduction)
    :returns: vector of virulence values through time
"""
function calculate_average_virulence()
    return sum(Np[:, 2:end] .* r1'; dims=2) ./ sum(Np[:, 2:end]; dims=2)
end

"""
Calculates the evenness of an Vector (shoutout to Pielou)
    :param n: vector
    :returns: evenness values through time
"""
function calculate_evenness(n)
    np = filter(x -> x > eps(), n)
    p = np ./ sum(np)
    ev = length(p) == 1 ? 0.0 : -sum(p .* log.(p)) * (1 / log(length(p)))
    return ev
end

# Plotting functions #########################################################
"""
Plots the numbers of individuals from each strain, comparing infected
and uninfected hosts as well as the total number of parasites.
    :param labels: vector of labels for each host
    :param plot_title: string of plot title
    :param upper_y_limit: upper limit for the y-axis
    :param upper_x_limit: upper limit for the x-axis
    :param png_path: path (as str) to write a PNG version of the plot
"""
function plot_population_numbers(
    labels::Array{String},
    plot_title::String,
    upper_y_limit::Int,
    upper_x_limit::Int,
    png_path::String
)
    # plot the number of infected hosts
    plot(
        Np[:, 2:end];
        c=:blue,
        lw=1.5,
        alpha=0.4,
        title=plot_title,
        xlabel="Time",
        ylabel="Number of individuals",
        label=labels,
        ylims=(0, upper_y_limit),
        xlims=(0, upper_x_limit),
    )
    # add the number of uninfected hosts
    plot!(
        Np[:, 1];
        c=:black,
        lw=1.5,
        label="Uninfected"
    )
    # add the total number of parasites
    plot!(
        sum(Np[:, 2:end]; dims=2);
        c=:red,
        label="Total parasites"
    )
    # save figure as a PNG
    png(png_path)
end

"""
Plots mortality data with regards to time.
    :param upper_y_limit: upper limit for the y-axis
    :param upper_x_limit: upper limit for the x-axis
    :param png_path: path (as str) to write a PNG version of the plot
"""
function plot_mortality(
    upper_y_limit::Int,
    upper_x_limit::Int,
    png_path::String
)
    plot(
        ui_w_avg;
        c=:black,
        lw=1.5,
        title="Average mortality in the population",
        xlabel="Time",
        ylabel="Mean mortality (ui)",
        leg=false,
        ylims=(0, upper_y_limit),
        xlims=(0, upper_x_limit)
    )
    # save figure as a PNG
    png(png_path)
end

"""
Plots the average reproductive rate (R0) with regards to time.
    :param upper_y_limit: upper limit for the y-axis
    :param upper_x_limit: upper limit for the x-axis
    :param png_path: path (as str) to write a PNG version of the plot
"""
function plot_reproductive_rate(
    upper_y_limit::Int,
    upper_x_limit::Int,
    png_path::String
)
    plot(
        H0_w + V0_w;
        c=:black,
        lw=1.5,
        title="Average R0 in the population",
        xlabel="Time",
        ylabel="Mean R0",
        leg=false,
        ylims=(0, upper_y_limit),
        xlims=(0, upper_x_limit),
    )
    # save figure as a PNG
    png(png_path)
end

"""
Plots the vertical basic reproductive ratio (V0) with regards to time.
    :param upper_y_limit: upper limit for the y-axis
    :param upper_x_limit: upper limit for the x-axis
    :param png_path: path (as str) to write a PNG version of the plot
"""
function plot_vertical_reproductive_ratio(
    upper_y_limit::Int,
    upper_x_limit::Int,
    png_path::String
)
    plot(
        V0_w;
        c=:black,
        lw=1.5,
        title="Average vertical cases in the population",
        xlabel="Time",
        ylabel="Mean V0",
        leg=false,
        ylims=(0, upper_y_limit),
        xlim=(0, upper_x_limit)
    )
    # save figure as a PNG
    png(png_path)
end

"""
Plots average horizontal transmission and virulence with regards to time.
    :param upper_y_limit: upper limit for the y-axis
    :param upper_x_limit: upper limit for the x-axis
    :param png_path: path (as str) to write a PNG version of the plot
"""
function plot_horizontal_and_virulence(
    upper_y_limit::Int,
    upper_x_limit::Int,
    png_path::String
)
    # plot the horizontal transmission
    plot(
        βy_w_avg;
        c=:black,
        title="Virulence and beta",
        label="Beta",
        xlabel="Time",
        ylabel="Mean virulence & \n Mean Beta",
        ylims=(0, upper_y_limit),
        xlims=(0, upper_x_limit),
    )
    # add virulence to the plot
    plot!(
        vir_w_avg;
        c=:blue,
        lw=1.5,
        label="Virulence"
    )
    # save figure as a PNG
    png(png_path)
end

"""
Calculates population evenness through its matrix.
Plots evenness data with regards to time.
    :param upper_y_limit: upper limit for the y-axis
    :param upper_x_limit: upper limit for the x-axis
    :param png_path: path (as str) to write a PNG version of the plot
"""
function plot_evenness(
    upper_y_limit::Int,
    upper_x_limit::Int,
    png_path::String
)
    # calculate the evenness through time
    evenness_data = mapslices(
        calculate_evenness,
        Np[:, 2:end];
        dims=2
    );

    # plots the evenness data
    plot(
        evenness_data,
        c=:black,
        lw=0.5,
        title="Evenness",
        xlabel="Time",
        ylabel="Relative abundance (log)",
        label=false,
        leg=false,
        ylims=(0, upper_y_limit),
        xlims=(0, upper_x_limit)
    )
    # save figure as a PNG
    png(png_path)
end


text:
function compute_dot(DX::Vector{Float64},DY::Vector{Float64})
assert(length(DX) == length(DY))
n = length(DX)
incx = incy = 1
product = ccall((:ddot, "libLAPACK"),
                 Float64,
                  (Ptr{Int32}, Ptr{Float64}, Ptr{Int32},
            Ptr{Float64}, Ptr{Int32}),
                  &n, DX, &incx, DY, &incy)
return product
end

max_stars_repo_path: Module 3/Chapter04/compute_dot.jl

max_issues_repo_path: Module 3/Chapter04/compute_dot.jl

max_forks_repo_path: Module 3/Chapter04/compute_dot.jl


text:
import POMDPs: MDP, Policy

#
# Define a policy that will act as an evolutionary/ genetic policy
#
mutable struct EvolutionPolicy <: Policy
    action_map::Vector{Int64}
end

#
# Default random intializer based on a games input
#
function EvolutionPolicy(pomdp::MDP; rng::AbstractRNG = MersenneTwister(0))
    default_map = rand(actions(pomdp), n_states(pomdp)) # all zeros
    return EvolutionPolicy(default_map)
end

#
# Default initializer with another policy as an input and a random
# mutate probability.
#
function EvolutionPolicy(pomdp::MDP, p1::EvolutionPolicy; p = 0.1, rng::AbstractRNG = MersenneTwister(0))
    p2 = EvolutionPolicy(pomdp; rng = rng)
    return EvolutionPolicy(pomdp, p1, p2; p = p, rng = rng)
end

#
# Default initializer that acts as a crossover between two different policies
#
function EvolutionPolicy(pomdp::MDP, p1::EvolutionPolicy, p2::EvolutionPolicy; p = 0.5, rng::AbstractRNG = MersenneTwister(0))

    l = n_states(pomdp)
    m = ifelse.(rand(rng, l) .< p, zeros(l), ones(l))

    new_action_map = m .* p1.action_map + (1 - m) .* p2.action_map

    return EvolutionPolicy(Vector{Int64}(new_action_map))
end

#
# Select and action from a policy based on a state
#
action(policy::EvolutionPolicy, s::Int64) = policy.action_map[s]

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
#初级传动比
function primGear_ratio(transm::PrimTransmission)
    drivenGear = transm.drivenGear  #离合器主动齿
    drivingGear = transm.drivingGear  #离合器从动齿
    return drivenGear/drivingGear  #初级传动比
end

#一档传动比
function firstGear_ratio(B::DrivenTransmission,A::DrivingTransmission)
    drivingGear = A.firstGear #主轴一档齿
    drivenGear = B.firstGear  #副轴一档齿
    return drivenGear/drivingGear  #一档传动比
end
#二档传动比
function secondGear_ratio(B::DrivenTransmission,A::DrivingTransmission)
    drivingGear = A.secondGear #主轴二档齿
    drivenGear = B.secondGear  #副轴二档齿
    return drivenGear/drivingGear  #二档传动比
end
#三档传动比
function thirdGear_ratio(B::DrivenTransmission,A::DrivingTransmission)
    drivingGear = A.thirdGear #主轴三档齿
    drivenGear = B.thirdGear  #副轴三档齿
    return drivenGear/drivingGear  #三档传动比
end
#四档传动比
function fourthGear_ratio(B::DrivenTransmission,A::DrivingTransmission)
    drivingGear = A.fourthGear #主轴四档齿
    drivenGear = B.fourthGear  #副轴四档齿
    return drivenGear/drivingGear  #四档传动比
end
#五档传动比
function fifthGear_ratio(B::DrivenTransmission,A::DrivingTransmission)
    drivingGear = A.fifthGear #主轴五档齿
    drivenGear = B.fifthGear  #副轴五档齿
    return drivenGear/drivingGear  #五档传动比
end
#六档传动比
function sixthGear_ratio(B::DrivenTransmission,A::DrivingTransmission)
    drivingGear = A.sixthGear #主轴六档齿
    drivenGear = B.sixthGear  #副轴六档齿
    return drivenGear/drivingGear  #六档传动比
end

#三级传动比计算
function wheel_ratio(wheel::Flywheel)
    return wheel.b_wheel/wheel.s_wheel

end

#发动机传动比
function vehiclespeed(enginename::String,ratio::Transmission; ca_speed = 7500, circumfeience = 1.97)
    ca_drvingGear = ratio.crankshaft.drivingGear
    ca_drvengGear = ratio.crankshaft.drivenGear

    drving_1_Gear = ratio.drivingshaft.firstGear
    drving_2_Gear = ratio.drivingshaft.secondGear
    drving_3_Gear = ratio.drivingshaft.thirdGear
    drving_4_Gear = ratio.drivingshaft.fourthGear
    drving_5_Gear = ratio.drivingshaft.fifthGear
    drving_6_Gear = ratio.drivingshaft.sixthGear

    drven_1_Gear = ratio.drivenshaft.firstGear
    drven_2_Gear = ratio.drivenshaft.secondGear
    drven_3_Gear = ratio.drivenshaft.thirdGear
    drven_4_Gear = ratio.drivenshaft.fourthGear
    drven_5_Gear = ratio.drivenshaft.fifthGear
    drven_6_Gear = ratio.drivenshaft.sixthGear

    primgearRatio =primGear_ratio(ratio.crankshaft)
    thirdgearRatio = wheel_ratio(ratio.flywheel)
    gear_1_ratio = firstGear_ratio(ratio.drivenshaft,ratio.drivingshaft)
    gear_2_ratio = secondGear_ratio(ratio.drivenshaft,ratio.drivingshaft)
    gear_3_ratio = thirdGear_ratio(ratio.drivenshaft,ratio.drivingshaft)
    gear_4_ratio = fourthGear_ratio(ratio.drivenshaft,ratio.drivingshaft)
    gear_5_ratio = fifthGear_ratio(ratio.drivenshaft,ratio.drivingshaft)
    gear_6_ratio = sixthGear_ratio(ratio.drivenshaft,ratio.drivingshaft)

    return VehicleSpeed(enginename,primgearRatio,[gear_1_ratio,gear_2_ratio,gear_3_ratio,gear_4_ratio,gear_5_ratio,gear_6_ratio],thirdgearRatio,circumfeience,ca_speed)

end

max_stars_repo_name: powerzf/VehicleCalculate

max_issues_repo_name: powerzf/VehicleCalculate

max_forks_repo_name: powerzf/VehicleCalculate


text:
# Autogenerated wrapper script for blis_jll for aarch64-linux-gnu
export blis

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("blis")
JLLWrappers.@declare_library_product(blis, "libblis.so.3")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        blis,
        "lib/libblis.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
##############   EXAMPLE-1    ##############
# Non-correlated errors in X and in Y

using LinearFitXYerrors
using Plots; gr()

# INPUT DATA:
# York, D. [1966] Least-squares fitting of a straight line. Canadian Journal of Physics, 44(5), pp.1079–1086
# Cantrell, C. [2008] Technical Note: Review of methods for linear least-squares fitting of data and application to atmospheric chemistry problems. Atmospheric Chem. & Physics, 8(17), pp.5477–5487

# Tables I and II from York [1966]

X = [0.0, 0.9, 1.8, 2.6, 3.3, 4.4, 5.2, 6.1, 6.5, 7.4]
Y = [5.9, 5.4, 4.4, 4.6, 3.5, 3.7, 2.8, 2.8, 2.4, 1.5]

# the weights ωX, ωY are the inverse of the variances
σX = 1 ./ sqrt.([1000., 1000, 500, 800, 200, 80,  60, 20, 1.8, 1])
σY = 1 ./ sqrt.([1., 1.8, 4, 8, 20, 20, 70, 70, 100, 500])
r = 0*X;

# COMPUTE and PLOT:
stxy = linearfitxy(X,Y; σX=σX, σY=σY, r=r, isplot=true);

# If assuming only erros in Y or in X:
sty = linearfitxy(X,Y; σX=0, σY=σY, r=r);
stx = linearfitxy(X,Y; σX=σX, σY=0, r=r);

dX = diff([extrema(X)...])[1]/7
x1, x2 = (-dX, dX) .+ extrema(X)
xx = [x1; X; x2]
plot!(xx, sty.a .+ sty.b*xx, color=:lime, lw=0.5, label="LinearFitXY (Y errors)")
plot!(xx, stx.a .+ stx.b*xx, color=:orange, lw=0.5, label="LinearFitXY (X errors)")

savefig("Example1_LinearFitXYerrors.png")
########################################


text:
using OMEinsum
using Test
using LinearAlgebra

@testset "OMEinsum.jl" begin
    include("Core.jl")
    include("EinRule.jl")
    include("utils.jl")
    include("einsum.jl")
    if Base.find_package("CuArrays") != nothing
        include("cueinsum.jl")
    end
    include("autodiff.jl")
    include("einsequence.jl")
    include("interfaces.jl")
end


text:
using GraffSDK
# Packages required for factor creation
using RoME
using Distributions

# 1. Import the initialization code.
cd(joinpath(Pkg.dir("GraffSDK"),"examples"))
include("0_Initialization.jl")

# 1a. Create a Configuration
config = loadGraffConfig("synchronyConfig.json")
#Create a hexagonal sessions
config.robotId = "DemoRobot"
config.sessionId = "DemoSession2"
println(getGraffConfig())

# --- Creating the robot if it doesn't exist ---- #
if !isRobotExisting()
    # Create a new one programatically - can also do this via the UI.
    println(" -- Robot '$(config.robotId)' doesn't exist, creating it...")
    newRobot = RobotRequest(config.robotId, "My New Bot", "Description of my neat robot", "Active")
    robot = addRobot(newRobot)
end

# --- Deleting the session if it exists---- #
if isSessionExisting()
    println(" -- Session '$(config.sessionId)' already exists for robot '$(config.robotId)', deleting and recreating it...")
    deleteSession()
end

# --- Make a new session - x0 will automatically be created during this call.

session = nothing
if isSessionExisting()
    println(" -- Session '$(config.sessionId)' already exists for robot '$(config.robotId)', retrieving it...")
    session = getSession()
else
    # Create a new one
    println(" -- Session '$(config.sessionId)' doesn't exist for robot '$(config.robotId)', creating it...")
    newSessionRequest = SessionDetailsRequest(config.sessionId, "A test dataset demonstrating basic graph creation.", "Pose2")
    session = addSession(newSessionRequest)
end

# --- Adding Variables ---- #
x1Request = VariableRequest("x1", "Pose2")
response = addVariable(x1Request)

x2Request = VariableRequest("x2", "Pose2")
response = addVariable(x2Request)

newLandmark = VariableRequest("l1", "Point2", ["LANDMARK"])
response = addVariable(newLandmark)

# --- Adding Factors ---- #

deltaMeasurement = [1.0; 1.0; pi/4]
pOdo = Float64[0.1 0 0; 0 0.1 0; 0 0 0.01]
p2p2 = Pose2Pose2(MvNormal(deltaMeasurement, pOdo.^2))
p2p2Conv = convert(PackedPose2Pose2, p2p2)
p2p2Request = FactorRequest(["x0", "x1"], "Pose2Pose2", p2p2Conv)

# Send the request for the x0->x1 link
addFactor(p2p2Request)
# Update the request to make the same link between x1 and x2
p2p2Request.variables = ["x1", "x2"]
addFactor(p2p2Request)

# Lastly, let's add the bearing+range factor between x1 and landmark l1
bearingDist = Normal(-pi/4, 0.1)
rangeDist = Normal(18, 1.0)
p2br2 = Pose2Point2BearingRange(bearingDist, rangeDist)
p2br2Conv = convert(PackedPose2Point2BearingRange, p2br2)
p2br2Request = FactorRequest(["x1", "l1"], "Pose2Point2BearingRange", p2br2Conv)
addFactor(p2br2Request)
# Now add the x1->l1 bearing+range factor
p2br2Request.variables = ["x2", "l1"]
addFactor(p2br2Request)

# We can now set these to ready and the server will solve it
putReady(true)


text:
mutable struct MyStruct a :: Int64 end
x = MyStruct(2)
x.a = true

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function send_stream(rd::IO, sock::TCPSocket, id::Int, mod)
    nb = bytesavailable(rd)
    if nb > 0
        d = read(rd, nb)
        s = String(copy(d))

        if !isempty(s)
            #info("sending $s to $(string(sock)) with id: $id")
            #= lock(serialize_lock) do
                remotecall_fetch(include_string, sock, Main, "
                    $(mod).print_to_console_remote(\"$(s)\", $(id))
                ")
            end =#
            lock(serialize_lock) do
                serialize(sock, (StdOutData(id,s),))
            end
        end
    end
end

function watch_stream(rd::IO, sock::TCPSocket, id::Int, mod)
    while !eof(rd) # blocks until something is available
        send_stream(rd, sock, id, mod)
        sleep(0.01) # a little delay to accumulate output
    end
end


text:
## General Market Model methods ##
function get_covariance(mm::AbstractMarketModel, i::Int)
  if isempty(mm.covariance)
    resize!(mm.covariance, mm.numberOfSteps)
    for j = 1:mm.numberOfSteps
      mm.covariance[j] = mm.pseudoRoots[j] * mm.pseudoRoots[j]'
    end
    i <= length(mm.covariance) || error("i must be less than or equal to covariance size")
  end

  return mm.covariance[i]
end

function total_covariance(mm::AbstractMarketModel, endIndex::Int)
  if isempty(mm.totalCovariance)
    resize!(mm.totalCovariance, mm.numberOfSteps)
    mm.totalCovariance[1] = get_covariance(mm, 1) # need to trigger calculation
    for j = 2:mm.numberOfSteps
      mm.totalCovariance[j] = mm.totalCovariance[j-1] + mm.covariance[j]
    end
    endIndex <= length(mm.covariance) || error("endIndex must be less than or equal to total covariance size")
  end

  return mm.totalCovariance[endIndex]
end


text:
@otestset "BVH" begin
    @testset "partition!" begin
        split = 0.5

        for i in 1:100000
            a = rand(5)
            b = copy(a)
            badresult::Bool = false

            lower, upper = partition!(a, (x) -> x, split)

            if lower !== nothing
                for i in lower
                    if i >= split
                        badresult = true
                        break
                    end
                end
            end

            if upper !== nothing
                for i in upper
                    if i <= split
                        badresult = true
                    end
                end
            end

            if badresult
                throw(ErrorException("array didn't partition: $(b)"))
            end

        end
    end
end # testset BVH