text:
# from S Das website, think adapting this should be the easiest thing to do
using SparseArrays;
using Serialization;
using LinearAlgebra;
using Arpack;
using KrylovKit;
using DelimitedFiles;
using RecursiveArrayTools;

include("Gaussian_kernel.jl");
include("Kernel_integral_Names.jl");

function Bistochastic_kernel(Data::Array{Float64,2}, offset::Int64, N::Int64, Q::Int64, epsilon::Float64, name::String="Expt", divide_by_Q::Bool=true  ) ::SparseMatrixCSC{Float64,Int64}
  K = Gaussian_kernel( Data, offset, N, Q, epsilon, divide_by_Q );
  serialize( open(Gaussian_kernel_name(name,N,Q,epsilon),"w"),K );
  Bistochastic_normalize!(K);
  serialize(open(Bistochastic_name(name,N,Q,offset,epsilon),"w"),K);
  #............Clear memory.................
  GC.gc();
  #.........................................
  return K;
end

function Bistochastic_normalize( K::SparseMatrixCSC{Float64,Int64}, DegreeVector_name::String = "DegVec", DensityVector_name::String = "DensVec") ::SparseMatrixCSC{Float64,Int64}
  N=size(K,1);
  #...............Compute d = K1_N.......................
  oneN=fill(1.0, N);
  d = K*oneN; d./=N; print("\nd computed");
  writedlm( DegreeVector_name, d, " ");
  #......................................................
  #............Compute q = K D^{-1} 1_N..................
  d_inv = d.^-1;
  D_inv = Diagonal(d_inv);
  q = K'*d_inv; q./=N; print("\nq computed");
  writedlm( DensityVector_name, q, " ");
  #......................................................
  #.............Bistochastic normalization...............
  q_sqrt_inv = q.^(-0.5);
  Qmatrix = Diagonal(q_sqrt_inv);
  tildeK = K*Qmatrix;
  tildeK.*= d_inv;
  tildeK./=sqrt(N); print("\nnormalization done");
  #......................................................
  q=0; d=0; d_inv=0; q_sqrt_inv=0; GC.gc();
  return tildeK;
end

"""
Stores ```1/m ∑_{j=1}^{m} ϕ_l(s_j)``` for ```1≤l≤L``` and stores in L-vector ***avg**. **KSX** is a ```m×N``` matrix storing ```KSX_{j,n} = k(s_j, x_n)``` for ```1≤j≤m, 1≤n≤N```. The ```{x_n ∈R^d}``` are the original data set, ```{s_j ∈R^d}``` are the new evaluation points. Matrix Γ= **Gamma_rescaled** is ```N×L```, with ```Γ_{n,l} = q(x_n)^{-0.5} γ_l(x_n) / λ_l```, where the ```γ_l``` are the right singular vectors in the bistochastic normalization.
"""
function Store_OS_avg_Bistoch( K_os::AbstractMatrix{Float64}, d_os :: Vector{Float64}, Gamma_rescaled::Matrix{Float64}, OS_avg :: Vector{Float64}, tv2 :: Vector{Float64}, tv3 :: Vector{Float64} )
	(m,N) = size(K_os);
	tv2 .= d_os.^-1; #tv2∈ R^m
	tv3.= K_os'*tv2; #tv3∈ R^N
	tv3./=(m*sqrt(N));
	OS_avg.= Gamma_rescaled'*tv3; #OS_avg∈ R^L
end

"""
**X** = ```[x_n 1≤n≤N] ∈R^{d×N}``` is the original data set, **S** = ```[s_j 1≤j≤m] ∈R^{d×m}``` are the evaluation points. **Grad** = ```J = [ J^{(j)} 1≤j≤m]^⊤``` is a ```dm×L``` matrix. This routine stores ```J^{(j)}_{:,l} = (∇ϕ_l)(s_j) ∈R^d```. **KSX** is a ```m×N``` matrix storing ```KSX_{j,n} = e(-|s_j-x_n|^2 / ϵ)``` for ```1≤j≤m, 1≤n≤N```. Matrix Γ= **Gamma_rescaled** is ```N×L```, with ```Γ_{n,l} = γ_l(x_n) / λ_l```, where the ```γ_l``` are the right singular vectors in the bistochastic normalization. **epsilon** = bandwidth parameter ϵ. **left** stores ```l(x_n)^{-1/2}```, produced in the bistochastic normalization.
"""
function Store_Gradient_Bistoch(X:: AbstractMatrix{Float64}, S::Matrix{Float64}, KSX::AbstractMatrix{Float64}, Gamma_rescaled::Matrix{Float64}, left::Vector{Float64}, epsilon:: Real, grad:: Matrix{Float64}, param:: Tuple )
  (d,m) = size(S);
  (N,L) = size(Gamma_rescaled);
  (tM1, tM2, tM3, tV1, tV2) = param;
  Store_line_sums(KSX, tV1, 2); tV1./=N; # tV1 \in R^{m} stores r(s_j)
  tM1.= Gamma_rescaled.*left; #tM1[n,l] = l(x_n)^{-1/2} γ_l(x_n) λ_l^{-1}
  grad.=0;
  for j = 1:m
    tM2.=-X;
    Threads.@threads for n in 1:N
      tM2[:,n] .+= S[:,j];
      tM2[:,n] .*= KSX[j,n];
    end
    #tM2[:,n] = exp( -|s_j - x_n|^2 / ϵ ) (s_j - x_n) ∈R^d
    Store_line_sums(tM2, tV2, 2); tV2./=N; # tV2 \in R^{d} stores (∇r)(s_j)
    tM3 .= tM2*tM1; #tM3 \in R^{d×L}
    tM3 ./= tV1[j];
    grad[ (j-1)*d+1:j*d, :] .= tM3;
    tM2 .= tV2*KSX[j,:]'; #tM3∈R^{d×N}
    tM3 .= tM2*tM1; #tM3 \in R^{d×L}
    tM3 ./= tV1[j]^2;
    grad[ (j-1)*d+1:j*d, :] .= tM3;
  end
  grad.*= 2/epsilon;
end

function RKHS_Koop_SVD( tildeK::SparseMatrixCSC{Float64,Int64}, L::Int64, name::String ) :: Array{Float64,2}
    #.............Constant parameters......................
    SVD_iterations = 10000;  SVD_tol=1E-5;
    #......................................................
    #.............Initializations..........................
    N=size(tildeK,1);#print("\nN = ",N," L = ",L);
    #......................................................
    #.....................SVD.............................
    Lambda,Phi,Gamma,info = KrylovKit.svdsolve( tildeK, L, krylovdim=2*L+1, verbosity=2, tol=SVD_tol, maxiter=SVD_iterations );
    L=size(Lambda,1);
    if L%2 == 0 L-=1; end
    print("\nL=",L);
    Phi=convert(Array,VectorOfArray(Phi[1:L]));
    print("\ntest : ", dot(Phi[:,1], Phi[:,4]));
    Gamma=convert(Array,VectorOfArray(Gamma[1:L]));
    writedlm(Gamma_name(name,N,Q,epsilon,L),Gamma," ");
    writedlm(Phi_name(name,N,Q,epsilon,L), Phi," ");
    writedlm(Lambda_name(name,N,Q,epsilon,L),Lambda[1:L]," ");
    writedlm(SVD_cnvrg_name(name,N,Q,epsilon,L),info.normres[1:L]," ");
    Threads.@threads for l in 1:L
      	Gamma[:,l] ./= Lambda[l];
    end
    writedlm(Gamma_rescaled_name(name,N,Q,epsilon,L),Gamma," ");
    print("\nnumber of conveged SVD pairs = ", info.converged, " after ",info.numiter," iterates, and ", info.numops," applications of the map");
    #display(info);
    #......................................................
    return Phi;
end

max_stars_repo_name: clairevalva/NLSA

max_issues_repo_name: clairevalva/NLSA

max_forks_repo_name: clairevalva/NLSA


text:
"""
    nsold01(x, FS, FPS, F!, J!=diffjac!; rtol=1.e-6, atol=1.e-12,
            maxit=20, dx=1.e-6, pdata=nothing)

This is Version .01. Nothing with a version number having a negative
exponent field can be trusted.

Nonlinear solvers from my books in Julia. This version has no globalization,
no quasi-Newton methods, and no Newton-Krylov. The mission here is to
duplicate the mixed precision results in my SIREV-ED submission.

# Inputs:
- x: initial iterate

- FS: Preallcoated storage for function. It is an N x 1 column vector

- FPS: preallcoated storage for Jacobian. It is an N x N matrix

- F!: function evaluation, the ! indicates that F! overwrites FS, your
    preallocated storage for the function.

- J!: Jacobian evaluation, the ! indicates that J! overwrites FPS, your
    preallocated storage for the Jacobian. If you leave this out the
    default is a finite difference Jacobian.

----------------------

# keyword arguments (kwargs):

- rtol and atol: relative and absolute error tolerances

- maxit: limit on nonlinear iterations

- dx: difference increment in finite-difference derivatives
      h=dx*norm(x)+1.e-8

- pdata: precomputed data for the function/Jacobian.
       Things will go better if you use this rather than hide the data
       in global variables within the module for your function/Jacobian

------------------------

# Using nsold.jl

Here are the rules as of June 6, 2019

F! is the nonlinear residual.
J! is the Jacobian evaluation.

I like to put everything in a module so that I can use global
(within the module) variables.

1) You allocate storage for the function and Jacobian in advance
   --> in the calling program <-- NOT in FS and FPS

FV=F!(FV,x) returns FV=F(x)

FP=J!(FV,FP,x) returns FP=F'(x);
    (FV,FP, x) must be the argument list, even if FP does not need FV.
    One reason for this is that the finite-difference Jacobian
    does and that is the default in the solver.

In the future J! will also be a matrix-vector product and FPS will
be the PREALLOCATED (!!) storage for the GMRES(m) Krylov vectors.

Lemme tell ya 'bout precision. I designed this code for full precision
functions and linear algebra in any precision you want. You can declare
FPS as Float64, Float32, or Float16 and nsold will do the right thing if
YOU do not destroy the declaration in your J! function. I'm amazed
that this works so easily.

If the Jacobian is reasonably well conditioned, I can see no reason
to do linear algebra in double precision

Don't try to evaluate function and Jacobian all at once because
that will cost you a extra function evaluation every time the line
search kicks in.

2) Any precomputed data for functions, Jacobians, matrix-vector products
   may live in global variables within a module containing F! and J!.
   Don't do that if you can avoid it. Use pdata instead.

"""
function nsold01(
    x,
    FS,
    FPS,
    F!,
    J! = diffjac!;
    rtol = 1.e-6,
    atol = 1.e-12,
    maxit = 20,
    dx = 1.e-6,
    pdata = nothing,
)

    EvalF!(FS, x, F!, pdata)
    resnorm = norm(FS)
    reshist = resnorm * ones(1,)
    tol = rtol * resnorm + atol
    itc = 0
    while resnorm > tol && itc < maxit
        #
        # FPF and FPS will now share the data for L and U, but are
        # different things. You can't say FPS=lu!(FPS) and then update
        # the Jacobian later in the nonlinear iteration.
        #
        # Don't believe me? Try it and see.
        #
        EvalJ!(FS, FPS, x, F!, J!, dx, pdata)
        FPF = lu!(FPS)
        step = -(FPF \ FS)
        x = x + step
        EvalF!(FS, x, F!, pdata)
        resnorm = norm(FS)
        push!(reshist, resnorm)
        itc = itc + 1
    end
    return (solution = x, ithist = reshist)
end

function EvalF!(FS, x, F!, pdata)
    F!
    if pdata == nothing
        F!(FS, x)
    else
        F!(FS, x, pdata)
    end
end

function EvalJ!(FS, FPS, x, F!, J!, dx, pdata)
    if J! != diffjac!
        if pdata == nothing
            J!(FS, FPS, x)
        else
            J!(FS, FPS, x, pdata)
        end
    else
        diffjac!(FS, FPS, F!, x, dx, pdata)
    end

end

function diffjac!(FS, FPS, F!, x, dx, pdata)
    precision = typeof(FPS[1, 1])
    h = dx * norm(x, Inf) + 1.e-8
    n = length(x)
    y = ones(size(x))
    FY = ones(size(x))
    for ic = 1:n
        y .= x
        y[ic] = y[ic] + h
        if pdata == nothing
            F!(FY, y)
        else
            F!(FY, y, pdata)
        end
        for ir = 1:n
            FPS[ir, ic] = (FY[ir] - FS[ir]) / h
        end
    end

end

ext: jl

max_stars_repo_name: ctkelley/MPResults


text:
using Test

@testset "Document" begin

    include("test_basic.jl")
    include("test_doc_changes.jl")
    include("test_utf16.jl")

end


text:
using Gnuplot, Random
x = randn(1000)
h = hist(x, bs = 0.1)
c = "#1E90FF"
@gp h.bins h.counts "with steps tit '' lc '$(c)' lw 2" "set grid"
save(term="pngcairo size 600,400", output = "stepPlot.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
#!/usr/bin/env julia

# Transverse field Ising model DMRG example.

using LatticeSweeper

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "-g"
        metavar = "g"
        help = "transverse field strength"
        arg_type = Float64
        required = true
    "-L"
        metavar = "L"
        help = "number of sites"
        arg_type = Int
        required = true
    "--max-sweeps"
        metavar = "S"
        help = "maximum number of sweeps"
        arg_type = Int
        required = true
    "--output-file"
        metavar = "PATH"
        help = "path to output file"
        arg_type = String
end
c = parse_args(ARGS, s, as_symbols=true)

g = c[:g]
L = c[:L]
max_sweeps = c[:max_sweeps]
output_file = c[:output_file]

# Pauli matrices.
id = [1.0 0.0; 0.0 1.0]
sigma_x = [0.0 1.0; 1.0 0.0]
sigma_z = [1.0 0.0; 0.0 -1.0]

println("[ ] Constructing Hamiltonian.")
H_tnsr = zeros(3, 3, 2, 2)
H_tnsr[1, 1, :, :] = id
H_tnsr[2, 1, :, :] = -g*sigma_x
H_tnsr[2, 2, :, :] = id
H_tnsr[2, 3, :, :] = -sigma_z
H_tnsr[3, 1, :, :] = sigma_z
@time H = MPO(H_tnsr, L)
println("[+] Constructed Hamiltonian.")

println("[ ] Constructing wavefunction.")
spin_up = [1.0, 0.0]
@time psi = MPS(spin_up, L)
println("[+] Constructed wavefunction.")

println("[ ] Sweeping.")
outputs = SweepOutput[]
if !isnothing(output_file)
    push!(outputs, SweepOutputFile(output_file))
end
push!(outputs, SweepOutputDynamic())
@time hist = dmrg!(psi, H, SweepSchedule(max_sweeps; tolerance=5e-12), outputs=outputs)
println("[+] Swept.")

hist.converged || @warn("Ground state not converged.")
if hist[end].mps_max_bond_dim >= hist[end].set.bond_max
    @warn("Bond dimension saturated at final sweep.")
end

# Ground state energy.
println_result("E0 = $(hist[end].energy)")
# Von Neumann entanglement entropy.
println_result("SvN = $(S_vn(hist[end].middle_eigvals))")


text:
using HPAT
using HPAT.API.LinearRegression
using MPI

@acc hpat function calcLinearRegression(file_name)
    points = DataSource(Matrix{Float64},HDF5,"/points", file_name)
    responses = DataSource(Vector{Float64},HDF5,"/responses", file_name)
    coeffs = LinearRegression(points, responses)
    return coeffs
end

c2 = calcLinearRegression(HPAT.getDefaultDataPath()*"/linear_regression_train_small.hdf5")
rank = MPI.Comm_rank(MPI.COMM_WORLD)
pes = MPI.Comm_size(MPI.COMM_WORLD)

MPI.Barrier(MPI.COMM_WORLD)
t1 = time_ns()
c2 = calcLinearRegression(HPAT.getDefaultDataPath()*"/linear_regression_train_large.hdf5")
t2 = time_ns()

if rank==0
	println("nodes: ", pes, "\nlinear regression exec time: ", (t2-t1)/1.0e9, "\nCoeffs: ", c2)
end

MPI.Barrier(MPI.COMM_WORLD)


text:
using LinearAlgebra

# Calculation of Intramolecular Force Fields from Second-Derivative Tensors
# Jorge M. Seminario

# Data struct
# Orca uses energy in Hartree units and length in Bohr
# 1 Hartee = 2625.5 kJ/mol; 1 Bohr = 0.0529177249 nm
# 1 [Hartee/Bohr^2] = 937582.935092 [kJ/mol/nm^2]
# 1 [Hartee/Bohr^2] = 2240.8815787 [kcal/mol/A^2]
#  xyz in angstrom; 1 A = 0.1 nm
struct Hess_data
    Hess::Matrix{Float64}
    xyz::Matrix{Float64}
    function Hess_data(Hess::Matrix{Float64}, xyz::Matrix{Float64})
        new(Hess*937582.935092 , 0.1*xyz)
    end
end

function bond_constant(atom1::Int64, atom2::Int64, MolSystem::Hess_data)
    # getting atom1-atom2 bond vector
    bond_vector = MolSystem.xyz[atom2, :] - MolSystem.xyz[atom1, :]
    bond_vector = bond_vector/norm(bond_vector)

    # getting submatrice of atom1-atom2 interaction
    HessAB = MolSystem.Hess[(3*atom1-2):(3*atom1), (3*atom2-2):3*atom2]
    eig, vec = eigen(-HessAB)

    # calculating bond force constant
    k = eig[1]*abs(dot(bond_vector, vec[:,1])) + eig[2]*abs(dot(bond_vector, vec[:,2])) + eig[3]*abs(dot(bond_vector, vec[:,3]))

    return(round(real(k), digits=3))

end

function angle_constant(atom1::Int64, atom2::Int64, atom3::Int64, MolSystem::Hess_data)
    # getting atom1-atom2; atom3-atom2 bond unit vectors ant their norms
    bond12_vector = MolSystem.xyz[atom2, :] - MolSystem.xyz[atom1, :]
    bond32_vector = MolSystem.xyz[atom2, :] - MolSystem.xyz[atom3, :]
    bond12_dist = norm(bond12_vector)
    bond32_dist = norm(bond32_vector)
    bond12_vector = bond12_vector/bond12_dist
    bond32_vector = bond32_vector/bond32_dist

    # getting submatrices of atom1-atom2 and atom2-atom3 interactions
    # calculating their eigen -values and -vectors
    Hess12 = MolSystem.Hess[(3*atom1-2):(3*atom1), (3*atom2-2):3*atom2]
    Hess32 = MolSystem.Hess[(3*atom3-2):(3*atom3), (3*atom2-2):3*atom2]
    eig12, vec12 = eigen(-Hess12)
    eig32, vec32 = eigen(-Hess32)

    # calculating direction vectors of small angle displacements
    normN_vector = cross(bond32_vector, bond12_vector)/norm(cross(bond32_vector, bond12_vector))
    disp12_vector = cross(normN_vector, bond12_vector)/norm(cross(normN_vector, bond12_vector))
    disp32_vector = cross(bond32_vector, normN_vector)/norm(cross(bond32_vector, normN_vector))

    # calculating angle force constant
    tempK1 = eig12[1]*abs(dot(disp12_vector, vec12[:,1])) + eig12[2]*abs(dot(disp12_vector, vec12[:,2])) + eig12[3]*abs(dot(disp12_vector, vec12[:,3]))
    tempK2 = eig32[1]*abs(dot(disp32_vector, vec32[:,1])) + eig32[2]*abs(dot(disp32_vector, vec32[:,2])) + eig32[3]*abs(dot(disp32_vector, vec32[:,3]))
    k = 1/(bond12_dist^2*tempK1) + 1/(bond32_dist^2*tempK2)
    k_theta = 1/k/2

    return(round(real(k_theta); digits=3))
end

function dihedral_constant(atom1::Int64, atom2::Int64, atom3::Int64, atom4::Int64,  MolSystem::Hess_data)
    # getting atom2-atom1; atom2-atom3; atom3-atom4; atom3-atom2 bond unit vectors ant their norms
    bond12_vector = MolSystem.xyz[atom2, :] - MolSystem.xyz[atom1, :]
    bond32_vector = MolSystem.xyz[atom2, :] - MolSystem.xyz[atom3, :]
    bond43_vector = MolSystem.xyz[atom3, :] - MolSystem.xyz[atom4, :]
    bond12_dist = norm(bond12_vector)
    bond32_dist = norm(bond32_vector)
    bond43_dist = norm(bond43_vector)
    bond12_vector = bond12_vector/bond12_dist
    bond32_vector = bond32_vector/bond32_dist
    bond43_vector = bond43_vector/bond43_dist
    bond23_vector = -bond32_vector/bond32_dist

    # getting submatrices of atom1-atom2 and atom3-atom4 interactions
    # calculating their eigen -values and -vectors
    Hess12 = MolSystem.Hess[(3*atom1-2):(3*atom1), (3*atom2-2):3*atom2]
    Hess34 = MolSystem.Hess[(3*atom3-2):(3*atom3), (3*atom4-2):3*atom4]
    eig12, vec12 = eigen(-Hess12)
    eig34, vec34 = eigen(-Hess34)

    # calculating unit normal vectors for planes atom1-atom2-atom3(ABC) and atom2-atom3-atom4(BCD)
    normNABC_vector = cross(bond32_vector, bond12_vector)
    normNBCD_vector = cross(bond43_vector, bond23_vector)
    normNABC_dist = norm(normNABC_vector)
    normNBCD_dist = norm(normNBCD_vector)
    normNABC_vector = normNABC_vector/normNABC_dist
    normNBCD_vector = normNBCD_vector/normNBCD_dist

    # calculating dihedral angle force constant
    tempK1 = eig12[1]*abs(dot(normNABC_vector, vec12[:,1])) + eig12[2]*abs(dot(normNABC_vector, vec12[:,2])) + eig12[3]*abs(dot(normNABC_vector, vec12[:,3]))
    tempK2 = eig34[1]*abs(dot(normNBCD_vector, vec34[:,1])) + eig34[2]*abs(dot(normNBCD_vector, vec34[:,2])) + eig34[3]*abs(dot(normNBCD_vector, vec34[:,3]))
    k = 1/((bond12_dist*normNABC_dist)^2*tempK1) + 1/((bond34_dist*normNBCD_dist)^2*tempK2)
    k_phi = 1/k

    return(abs(real(k_phi)))
end

function improper_constant(atom1::Int64, atom2::Int64, atom3::Int64, atom4::Int64,  MolSystem::Hess_data)
    # getting atom1-atom2; atom1-atom3; atom1-atom4 and atom4-atom3; atom2-atom3 bond unit vectors ant their norms
    # getting coordinates of central atom1 and atom2
    Atom1_point = MolSystem.xyz[atom1, :]
    Atom2_point = MolSystem.xyz[atom2, :]
    bond12_vector = Atom2_point - Atom1_point
    bond13_vector = MolSystem.xyz[atom3, :] - Atom1_point
    bond14_vector = MolSystem.xyz[atom4, :] - Atom1_point
    bond43_vector = MolSystem.xyz[atom3, :] - MolSystem.xyz[atom4, :]
    bond23_vector = MolSystem.xyz[atom3, :] - Atom2_point
    bond12_dist = norm(bond12_vector)
    bond13_dist = norm(bond13_vector)
    bond14_dist = norm(bond14_vector)
    bond43_dist = norm(bond43_vector)
    bond23_dist = norm(bond23_vector)
    bond12_vector = bond12_vector/bond12_dist
    bond13_vector = bond13_vector/bond13_dist
    bond14_vector = bond14_vector/bond14_dist
    bond43_vector = bond43_vector/bond43_dist
    bond23_vector = bond23_vector/bond23_dist

    # getting submatrices of atom1-atom2 and atom1-atom3 and atom1-atom4 interactions
    # calculating their eigen -values and -vectors
    Hess12 = MolSystem.Hess[(3*atom1-2):(3*atom1), (3*atom2-2):3*atom2]
    Hess13 = MolSystem.Hess[(3*atom1-2):(3*atom1), (3*atom3-2):3*atom3]
    Hess14 = MolSystem.Hess[(3*atom1-2):(3*atom1), (3*atom4-2):3*atom4]
    eig12, vec12 = eigen(-Hess12)
    eig13, vec13 = eigen(-Hess13)
    eig14, vec14 = eigen(-Hess14)

    # calculating unit normal vectors for planes atom2-atom3-atom4(BCD)
    normNBCD_vector = cross(bond43_vector, bond23_vector)
    normNBCD_dist = norm(normNBCD_vector)
    normNBCD_vector = normNBCD_vector/normNBCD_dist

    # calculating improper dihedral angle force constant

    # calculating in
    tempK1 = eig12[1]*abs(dot(normNBCD_vector , vec12[:,1])) + eig12[2]*abs(dot(normNBCD_vector , vec12[:,2])) + eig12[3]*abs(dot(normNBCD_vector , vec12[:,3]))
    tempK2 = eig13[1]*abs(dot(normNBCD_vector , vec13[:,1]))+ eig13[2]*abs(dot(normNBCD_vector , vec13[:,2])) + eig13[3]*abs(dot(normNBCD_vector , vec13[:,3]))
    tempK3 = eig14[1]*abs(dot(normNBCD_vector , vec14[:,1])) + eig14[2]*abs(dot(normNBCD_vector , vec14[:,2])) + eig14[3]*abs(dot(normNBCD_vector , vec14[:,3]))

    # calculating projection of altitude of triangle ABC into plane BCD
    # determinating altitude vector of altitude through system of linear equations
    # Statement 1: equation of line BC; Statement 2: dot product is 0 for (BC*altitude),  MatrCoef*r = x
    # BC is bond23_vector; Atom2_point is B point; Atom1_point is A point

    MatrCoef = [1/bond23_vector[1]   -1/bond23_vector[2]           0;
                          0           1/bond23_vector[2]   -1/bond23_vector[3];
                bond23_vector[1]        bond23_vector[2]    bond23_vector[3]]

    x = [(Atom2_point[1]/bond23_vector[1] - Atom2_point[2]/bond23_vector[2]);
         (Atom2_point[2]/bond23_vector[2] - Atom2_point[3]/bond23_vector[3]);
        -(dot(Atom1_point, bond23_vector))]

    HABC_vector = Atom1_point - MatrCoef\x
    normHABCD_dist = norm(HABC_vector - HABC_vector*dot(HABC_vector, normNBCD_vector))
    k_omega = normHABCD_dist^2*(tempK1 + tempK2 + tempK3)

    return(abs(real(k_omega)))
end

max_stars_repo_name: Dimitrifulgur/SeminarioFUERZA

max_issues_repo_name: Dimitrifulgur/SeminarioFUERZA

max_forks_repo_name: Dimitrifulgur/SeminarioFUERZA


text:
# Base for Arm models

abstract type BanditArmBase end

function pull!( arm::BanditArmBase )
    error( "No Implementation of pull!() for ", typeof(arm) )
end

function tick!( arm::BanditArmBase )
    error( "No Implementation of tick!() for ", typeof(arm) )
end

function reset!( arm::BanditArmBase )
    error( "No Implementation of reset!() for ", typeof(arm) )
end


text:
using Chemostat_Nanchen2006
using Test

@testset "Chemostat_Nanchen2006.jl" begin
    # Write your tests here.
end


text:
using WassersteinDictionaries
using Documenter

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

makedocs(;
    modules=[WassersteinDictionaries],
    authors="Tobias M. Blickhan, Michael Kraus",
    repo="https://github.com/JuliaRCM/WassersteinDictionaries.jl/blob/{commit}{path}#{line}",
    sitename="WassersteinDictionaries.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaRCM.github.io/WassersteinDictionaries.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaRCM/WassersteinDictionaries.jl",
    devbranch="main",
)

max_stars_repo_name: JuliaRCM/WassersteinDictionaries

max_issues_repo_name: JuliaRCM/WassersteinDictionaries

max_forks_repo_name: JuliaRCM/WassersteinDictionaries


text:
import AQFED.Math: norminv
using Statistics
import AQFED.TermStructure: HestonModel, discountFactor, forward
import AQFED.Random: next!, nextn!, skipTo

function ndims(model::HestonModel, specificTimes::Vector{Float64}, timestepSize::Float64)
    genTimes = pathgenTimes(model, specificTimes, timestepSize)
    return (length(genTimes) - 1) * 2 #0.0 does not count
end

#DVSS2X discretization scheme for Heston
function simulateDVSS2X(
    rng,
    model::HestonModel{T},
    spot::Float64,
    payoff::VanillaOption,
    start::Int,
    nSim::Int,
    timestepSize::Float64;
    withBB = false,
    cacheSize = 0
) where {T}
    specTimes = specificTimes(payoff)
    tte = specTimes[end]
    df = discountFactor(model, tte)
    genTimes = pathgenTimes(model, specTimes, timestepSize)
    logpathValues = Vector{T}(undef, nSim)
    local bb, cache
    if withBB
        bb = BrownianBridgeConstruction(genTimes[2:end])
        if cacheSize == 0
            nSteps = length(genTimes) - 1
            cacheSize = ceil(Int, log2(nSteps)) * 4
        end
        # cache = BBCache{Int,Array{Float64,2}}(cacheSize)
        cache = BBCache{Int,Vector{Float64}}(cacheSize)
    end

    u1 = Vector{Float64}(undef, nSim)
    u2 = Vector{Float64}(undef, nSim)
    xsi = u1
    v = Vector{T}(undef, nSim)
    t0 = genTimes[1]
    lnspot = log(spot)
    lnf0 = logForward(model, lnspot, t0)
    logpathValues .= lnf0
    v .= model.v0
    ρBar = sqrt(1 - model.ρ^2)
    ndimsh = length(genTimes) - 1
    local payoffValues
    for (dim, t1) in enumerate(genTimes[2:end])
        h = t1 - t0
        ekdth = exp(-model.κ * h / 2)
        c0 = model.θ * h / 4
        c1 = (ekdth - 1) / (model.κ * 2)
        c2 = model.θ * (1 - ekdth)
        if withBB
            d = dim
            transformByDim!(bb, rng, start, d, xsi, cache)
            skipTo(rng, d + ndimsh, start)
            next!(rng, d + ndimsh, u2)
        else
            d = 2 * dim - 1
            skipTo(rng, d, start)
            nextn!(rng, d, xsi)
            xsi .*= sqrt(h)
            skipTo(rng, d + 1, start)
            next!(rng, d + 1, u2)
        end

        vTmp = @. (v * ekdth + c2) / model.σ^2
        yr = map(
            (vTmp, u2) -> ifelse(
                vTmp > 1.896 * h,  #limit for y1 > 0 => 1.89564392373896
                dvss2_case1(vTmp, u2, h),
                dvss2_case2(vTmp, u2, h),
            ),
            vTmp,
            u2,
        )
        xb = vTmp # reuse var
        @. xb =
            logpathValues - c0 +
            c1 * (v - model.θ) +
            model.σ * (ρBar * xsi * sqrt((vTmp + yr) / 2) + model.ρ * (yr - vTmp))
        yr .= model.σ^2 * yr
        lnf1 = logForward(model, lnspot, t1)
        @. logpathValues = lnf1 - lnf0 + xb - c0 + c1 * (yr - model.θ)
        @. v = yr * ekdth + c2
        if t1 == tte
            pathValues = xb #reuse Var
            @. pathValues = exp(logpathValues)
            payoffValues = map(x -> evaluatePayoff(payoff, x, df), pathValues)
        end
        t0 = t1
        lnf0 = lnf1
    end
    payoffMean = mean(payoffValues)
    return payoffMean, stdm(payoffValues, payoffMean) / sqrt(length(payoffValues))
end

@inline function dvss2_case1(vTmp, u2, h)
    delta = (21.0 / 4 * h + 12 * vTmp) * h
    sqrtDelta = sqrt(delta)
    s = 1.5 * h
    twovh = 2 * vTmp * h / sqrtDelta
    p1 = twovh / (sqrtDelta - s)
    p2 = twovh / (sqrtDelta + s)
    y0 = vTmp
    y1 = vTmp + (s - sqrtDelta) / 2
    y2 = y1 + sqrtDelta
    if (u2 < p1)
        return y1
    elseif u2 > 1 - p2
        return y2
    else
        return y0
    end
end

@inline function dvss2_case2(vTmp, u2, h)
    denom = (2 * vTmp + h)
    s = denom + (5 * vTmp + 2 * h) * h / denom
    delta = h * (((16 * vTmp + 33 * h) * vTmp + 18 * h^2) * vTmp + 3 * h^3)
    sqrtDelta = sqrt(delta) / denom
    p1 = vTmp * (2 * (vTmp + h) - s - sqrtDelta) / (sqrtDelta * (sqrtDelta - s))
    p2 = vTmp * (2 * (vTmp + h) - s + sqrtDelta) / (sqrtDelta * (sqrtDelta + s))
    y0 = 0.0
    y1 = (s - sqrtDelta) * 0.5
    y2 = y1 + sqrtDelta
    if (u2 < p1)
        return y1
    elseif u2 > 1 - p2
        return y2
    else
        return y0
    end
end

#Euler Full Truncation discretization scheme for Heston
function simulateFullTruncation(
    rng,
    model::HestonModel{T},
    spot::Float64,
    payoff::VanillaOption,
    start::Int,
    nSim::Int,
    timestepSize::Float64;
    withBB = false,
    cacheSize = 0
) where {T}
    specTimes = specificTimes(payoff)
    tte = specTimes[end]
    df = discountFactor(model, tte)
    genTimes = pathgenTimes(model, specTimes, timestepSize)
    logpathValues = Vector{T}(undef, nSim)
    local bb, cache
    if withBB
        bb = BrownianBridgeConstruction(genTimes[2:end])
        if cacheSize == 0
            nSteps = length(genTimes) - 1
            cacheSize = ceil(Int, log2(nSteps)) * 4
        end
        cache = BBCache{Int,Array{Float64,2}}(cacheSize)
        # cache = BBCache{Int,Vector{Float64}}(cacheSize)
    end

    u = Array{Float64}(undef, (nSim, 2))
    u1 = @view u[:, 1]
    u2 = @view u[:, 2]
    v = Vector{T}(undef, nSim)
    sqrtmv = Vector{T}(undef, nSim)
    t0 = genTimes[1]
    lnspot = log(spot)
    lnf0 = logForward(model, lnspot, t0)
    logpathValues .= lnf0
    v .= model.v0
    ρBar = sqrt(1 - model.ρ^2)
    ndimsh = length(genTimes) - 1
    local payoffValues
    for (dim, t1) in enumerate(genTimes[2:end])
        h = t1 - t0
        sqrth = sqrt(h)
        if withBB
            transformByDim!(bb, rng, start, dim, u, cache)
        else
            d = 2 * dim - 1
            skipTo(rng, d, start)
            nextn!(rng, d, u1)
            skipTo(rng, d + 1, start)
            nextn!(rng, d + 1, u2)
            #@inbounds @. u *= sqrth
            u1 .*= sqrth
            u2 .*= sqrth
        end
        lnf1 = logForward(model, lnspot, t1)
        @. sqrtmv = sqrt(max(v, 0))
        @. logpathValues +=
            lnf1 - lnf0 - (0.5 * sqrtmv^2) * h + sqrtmv * (u1 * ρBar + u2 * model.ρ)
        # for (i,p) in enumerate(logpathValues)
        #     if isnan(p)
        #         println(i," nan ",sqrtmv[i], " ",u1[i]," ", u2[i])
        #     end
        # end
        @. v += model.κ * (model.θ - sqrtmv^2) * h + model.σ * sqrtmv * u2

        if t1 == tte
            pathValues = sqrtmv #reuse Var
            @. pathValues = exp(logpathValues)
            payoffValues = map(x -> evaluatePayoff(payoff, x, df), pathValues)
        end
        t0 = t1
        lnf0 = lnf1
    end
    payoffMean = mean(payoffValues)
    return payoffMean, stdm(payoffValues, payoffMean) / sqrt(length(payoffValues))
end

function simulateFullTruncationIter(
    rng,
    model::HestonModel{T},
    spot::Float64,
    payoff::VanillaOption,
    start::Int,
    nSim::Int,
    timestepSize::Float64;
    withBB = false
) where {T}
    specTimes = specificTimes(payoff)
    tte = specTimes[end]
    df = exp(-model.r * tte)
    genTimes = pathgenTimes(model, specTimes, timestepSize)
    local bb
    if withBB
        bb = BrownianBridgeConstruction(genTimes[2:end])
    end
    lnspot = log(spot)
    lnf0 = logForward(model, lnspot, tte)
    ρBar = sqrt(1 - model.ρ^2)
    z = Vector{Float64}(undef, (length(genTimes) - 1) * 2)
    u = Array{Float64}(undef, (2, length(genTimes) - 1))
    skipTo(rng, start)
    payoffMean = 0.0
    for sim = 1:nSim
        if withBB
            nextn!(rng, z)
            transform!(bb, z, u)
        else
            nextn!(rng, z)
            t0 = genTimes[1]
            @inbounds for i = 1:length(genTimes)-1
                dim = 2 * i - 1
                t1 = genTimes[i+1]
                u[1, i] = z[dim] * sqrt(t1 - t0)
                u[2, i] = z[dim+1] * sqrt(t1 - t0)
                t0 = t1
            end
        end
        t0 = genTimes[1]
        logpathValue = lnf0
        v = model.v0
        local payoffValue
        @inbounds for i = 1:length(genTimes)-1
            t1 = genTimes[i+1]
            # dim = 2 * i - 1
            u1 = u[1, i]
            u2 = u[2, i]
            h = t1 - t0
            sqrtmv = sqrt(max(v, 0))
            lnf1 = logForward(model, lnspot, t1)
            logpathValue += lnf1 - lnf0 - 0.5 * sqrtmv^2 * h + sqrtmv * (u1 * ρBar + u2 * model.ρ)
            v += model.κ * (model.θ - sqrtmv^2) * h + model.σ * sqrtmv * u2

            if t1 == tte
                pathValue = exp(logpathValue)
                payoffValue = evaluatePayoff(payoff, pathValue, df)
            end
            t0 = t1
            lnf0 = lnf1
        end

        payoffMean += payoffValue
    end
    payoffMean /= nSim
    return payoffMean, NaN
end


text:
  # differential evolution strategy
  using Distributions
  using LinearAlgebra

  function des(par, fn, lower, upper, control)
    @time begin
     function controlParam(name, default)
       if haskey(control, name) == false
           return default
       else
           return control[name]
       end
     end

     function drop(M, r=nothing, c=nothing)
        dr = collect(1:size(M,1))
        dc = collect(1:size(M,2))
        isnothing(r) ? nothing : splice!(dr,r)
        isnothing(c) ? nothing : splice!(dc,c)
        M[dr,dc]
      end

     function sampleFromHistory(history, historySample, lambda)
       ret = []
       for i in 1:lambda
         x = sample(1:size(history[convert(Int64,historySample[i])],2), 1)
         append!(ret, x)
       end
       return ret
     end

     function deleteInfsNaNs(x)
       for i in eachindex(x)
         if isnan(x[i])
           x[i] = prevfloat(typemax(Float64))
         end
       end
       for i in eachindex(x)
         if isinf(x[i])
           x[i] = prevfloat(typemax(Float64))
         end
       end
       return x
     end

     function fn_(x)
       if all(x.>=lower) && all(x.<=upper)
         counteval = counteval + 1
         return fn(x)
       else
         return prevfloat(typemax(Float64))
       end
     end

    function apply(x, f)
      r = []
      for i in 1:size(x, 2)
        append!(r, f(x[:,i]))
      end
      return r
    end

   function fn_l(P)
     if ndims(P) == 2 && size(P, 2) != 1 #if is.matrix()
        if (counteval + size(P, 2)) <= budget
          x = apply(P, fn_)
          return x
        else
          ret = []
          budLeft = budget - counteval
          if budLeft > 0
            for i in 1:budLeft
              append!(ret, fn_(P[:,i]))
            end
          end
          x = append!(ret, repeat([prevfloat(typemax(Float64))], size(P, 2) - budLeft))
          return x
        end
     else
       if counteval < budget
         return fn_(P)
       else
         return prevfloat(typemax(Float64))
       end
     end
   end

   function fn_d(P, P_repaired, fitness)
     P = deleteInfsNaNs(P)
     P_repaired = deleteInfsNaNs(P_repaired)

     if ndims(P) == 2 && ndims(P_repaired) == 2
       repairedInd = []
       for i in 1:size(P,2)
         append!(repairedInd,  all(P[:,i].!= P_repaired[:,i]))
       end
       P_fit = fitness
       P_pom = (P - P_repaired).^2
       vecDist = apply(P_pom,sum)
       if length(vecDist) == 1
         P_fit[findall(repairedInd)] = vecDist[findall(repairedInd)] + worst_fit
       else
         P_fit[findall(repairedInd)] = vecDist[findall(repairedInd)].+ worst_fit
       end
       P_fit = deleteInfsNaNs(P_fit)
       return P_fit
     else
       P_fit = fitness
       if P != P_repaired
        P_pom = (P - P_repaired).^2
         P_fit = worst_fit + sum(P_pom)
         P_fit = deleteInfsNaNs(P_fit)
       end
       return P_fit
     end
   end

   function bounceBackBoundary2(x)
     pom = true
     for i in eachindex(x)
       if x[i] <= lower || x[i] >= upper
         pom = false
       end
     end
     if pom == true
       return x
     elseif any(x.<lower)
       for i in findall(x.<lower)
         x[i] = lower + abs(lower - x[i]) % (upper - lower)
       end
     elseif any(x.>upper)
       for i in findall(x.>upper)
         x[i] = upper - abs(upper - x[i]) % (upper - lower)
       end
     end
     x = deleteInfsNaNs(x)
     return bounceBackBoundary2(x)
   end

   N = length(par)

   # Algorithm parameters
   Ft          = controlParam("Ft", 1) # Scaling factor of difference vectors (a variable!)
   initFt      = controlParam("initFt", 1)
   stopfitness = controlParam("stopfitness", -Inf) # Fitness value after which the convergence is reached
   # Strategy parameter setting:
   budget      = controlParam("budget", 10000 * N) # The maximum number of fitness function calls
   initlambda  = controlParam("lambda", 4 * N) # Population starting size
   lambda      = initlambda ## Population size
   mu          = controlParam("mu", floor(lambda / 2)) # Selection size
   weights     = controlParam("weights", Base.log.(repeat([mu + 1],convert(Int64, mu))) - Base.log.(1:convert(Int64, mu))) # Weights to calculate mean from selected individuals
   weights     = weights / sum(weights) # weights are normalized by the sum
   weightsSumS = sum(weights.^2) # weights sum square
   mueff       = controlParam("mueff", sum(weights)^2 / sum(weights.^2)) # Variance effectiveness factor
   cc          = controlParam("ccum", mu / (mu + 2)) # Evolution Path decay factor
   pathLength  = controlParam("pathLength", 6) # Size of evolution path
   cp          = controlParam("cp", 1 / sqrt(N)) # Evolution Path decay factor
   maxiter     = controlParam("maxit", floor(budget / (lambda + 1))) # Maximum number of iterations after which algorithm stops
   c_Ft        = controlParam("c_Ft", 0)
   pathRatio   = controlParam("pathRatio", sqrt(pathLength)) # Path Length Control reference value
   histSize    = controlParam("history", ceil(6 + ceil(3 * sqrt(N)))) # Size of the window of history - the step length history
   Ft_scale    = controlParam("Ft_scale", ((mueff + 2) / (N + mueff + 3)) / (1 + 2 * max(0, sqrt((mueff - 1) / (N + 1)) - 1) + (mueff + 2) / (N + mueff + 3)))
   tol         = controlParam("tol", 10^-12)
   counteval   = 0 # Number of function evaluations
   sqrt_N      = sqrt(N)

   log_all = controlParam("diag", false)
   log_Ft = controlParam("diag.Ft", log_all)
   log_value = controlParam("diag.value", log_all)
   log_mean = controlParam("diag.mean", log_all)
   log_meanCord = controlParam("diag.meanCords", log_all)
   log_pop = controlParam("diag.pop", log_all)
   log_bestVal = controlParam("diag.bestVal", log_all)
   log_worstVal = controlParam("diag.worstVal", log_all)
   log_eigen = controlParam("diag.eigen", log_all)
   Lamarckism = controlParam("Lamarckism", false)

    best_fit = Inf # The best fitness found so far
    best_par = nothing # The best solution found so far
    worst_fit = nothing # The worst solution found so far
    last_restart = 0
    restart_length = 0
    restart_number = 0

    if log_Ft
      Ft_log = zeros(0, 1)
    end
    if log_value
      value_log = zeros(0, lambda)
    end
    if log_mean
      mean_log = zeros(0, 1)
    end
    if log_meanCord
      meanCords_log = zeros(0, N)
    end
    if log_pop
      pop_log = zeros(N, lambda, maxiter)
    end
    if log_bestVal
      bestVal_log = zeros(0, 1)
    end
    if log_worstVal
      worstVal_log = zeros(0, 1)
    end
    if log_eigen
      eigen_log = zeros(0, N)
    end

    # buffers:
    dMean = zeros(N, convert(Int64,histSize))
    FtHistory = zeros(convert(Int64,histSize)) # Array buffer containing 'histSize' last values of 'Ft'
    pc = zeros(N, convert(Int64,histSize))

    # Initialize internal strategy parameters
    msg = nothing # Reason for terminating
    restart_number = -1
    iter = 0 # Number of iterations

    while counteval < budget

      restart_number = restart_number + 1
      mu = floor(lambda / 2)
      weights = Base.log.(repeat([mu + 1],convert(Int64, mu))) - Base.log.(1:convert(Int64, mu))
      weights = weights / sum(weights)
      weightsPop = Base.log.(repeat([lambda + 1],convert(Int64, lambda))) - Base.log.(1:convert(Int64, lambda))
      weightsPop = weightsPop / sum(weightsPop)

      histHead = 0 # Pointer to the history buffer head
      iter = 0 # Number of iterations
      history = [] # List stores best 'mu'(variable) individuals for 'hsize' recent iterations
      Ft = initFt

      # Create fisrt population
      population = rand(N,lambda).*(0.8*upper-0.8*lower).+0.8*lower

      cumMean = (upper + lower) / 2

      populationRepaired = reshape(apply(population, bounceBackBoundary2), N, lambda)

      if Lamarckism == true
        population = populationRepaired
      end

      selection = fill(0, convert(Int64,mu))
      selectedPoints =  zeros(N, convert(Int64,mu))
      fitness = fn_l(population)

      oldMean = zeros(N)
      newMean = par
      limit = 0
      worst_fit = maximum(fitness)

      # Store population and selection means
     popMean = drop(population * weightsPop)
     muMean = newMean

     # Matrices for creating diffs
     diffs = zeros(N, lambda)
     x1sample = zeros(lambda)
     x2sample = zeros(lambda)

     chiN = sqrt(N)

     histNorm = 1 / sqrt(2)
     counterRepaired = 0

     stoptol = false

      while counteval < budget && !stoptol
        iter = iter + 1
        histHead = (histHead % histSize) + 1

        mu = floor(lambda / 2)
        weights = Base.log.(repeat([mu + 1],convert(Int64, mu))) - Base.log.(1:convert(Int64, mu))
        weights = weights / sum(weights)

        if log_Ft
           Ft_log = vcat(Ft_log, Ft)
        end
        if log_value
           value_log = vcat(value_log, fitness)
        end
        if log_mean
          mean_log = vcat(mean_log, fn_l(bounceBackBoundary2(newMean)))
        end
        if log_meanCord
          meanCords_log = vcat(meanCords_log, newMean)
        end
        if log_pop
          pop_log[:, :, iter] = population
        end
        if log_bestVal
           bestVal_log = vcat(bestVal_log, minimum(minimum(bestVal_log), minimum(fitness)))
        end
        if log_worstVal
           worstVal_log = vcat(worstVal_log, maximum(maximum(worstVal_log), maximum(fitness)))
        end
        if log_eigen
           eigen_log = vcat(eigen_log, reverse(sort(eigvals(cov(transpose(population))))))
        end

       # Select best 'mu' individuals of popu-lation
       selection = sortperm(vec(fitness))[1:convert(Int64,mu)]
       selectedPoints = population[:, selection]

       push!(history, selectedPoints.* histNorm./ Ft)

       # Calculate weighted mean of selected points
       oldMean = newMean
       newMean = drop(selectedPoints * weights)

       # Write to buffers
       muMean = newMean
       dMean[:, convert(Int64,histHead)] = (muMean - popMean)./ Ft

       step = (newMean - oldMean)./ Ft

       # Update Ft
       FtHistory[convert(Int64,histHead)] = Ft
       oldFt = Ft

      # Update parameters
      if  convert(Int64,histHead) == 1
        pc[:, convert(Int64,histHead)] = fill(0.0, N).*(1 - cp)./ sqrt(N) + step.*sqrt(mu * cp * (2 - cp))
      else
        pc[:, convert(Int64,histHead)] = pc[:, convert(Int64,histHead - 1)].*(1 - cp) + step.*sqrt(mu * cp * (2 - cp))
      end

      # Sample from history with uniform distribution
      if iter < histSize
        limit = histHead
      else
        limit = histSize
      end
      historySample = sample(1:limit, lambda, replace=true)
      historySample2 = sample(1:limit, lambda, replace=true)

      x1sample = sampleFromHistory(history, historySample, lambda)
      x2sample = sampleFromHistory(history, historySample, lambda)

      # Make diffs
      for i in 1:lambda
        x1 = history[convert(Int64,historySample[i])][:, x1sample[i]]
        x2 = history[convert(Int64,historySample[i])][:, x2sample[i]]

        diffs[:, i].= sqrt(cc) * ((x1 - x2) + dMean[:, convert(Int64,historySample[i])]).*randn(1) + pc[:, convert(Int64,historySample2[i])].*sqrt(1 - cc).* randn(1)[1]
      end

     # New population
     population = newMean.+ diffs.*Ft.+ reshape(tol * (1 - 2 / N^2)^(iter / 2) * rand(Normal(), length(diffs)), N, lambda)./ chiN
     population = deleteInfsNaNs(population)

     # Check constraints violations
     # Repair the individual if necessary
     populationTemp = population
     populationRepaired = apply(population, bounceBackBoundary2) #mapslices(bounceBackBoundary2, population, dims=[1])
     populationRepaired = reshape( hcat(populationRepaired...), size(populationTemp,1), size(populationTemp, 2))
     counterRepaired = 0

     for tt in 1:size(populationTemp,2)
       if any(populationTemp[:, tt] != populationRepaired[:, tt])
         counterRepaired = counterRepaired + 1
       end
     end

     if Lamarckism == true
        population = populationRepaired
     end

     popMean =  drop(population * weightsPop)

      # Evaluation
     fitness = fn_l(population)

     if Lamarckism == false
        fitnessNonLamarcian = fn_d(population, populationRepaired, fitness)
     end

     wb = argmin(vec(fitness))

     if fitness[wb] < best_fit[1]
        best_fit = fitness[wb]
        if Lamarckism == true
          best_par = population[:, wb]
        else
          best_par = populationRepaired[:, wb]
        end
     end

      # Check worst fit:
     ww = argmax(fitness)
     if fitness[ww] > worst_fit
        worst_fit = fitness[ww]
     end

      # Fitness with penalty for nonLamarcian approach
      if Lamarckism == false
        fitness = fitnessNonLamarcian
      end

      # Check if the middle point is the best found so far
      cumMean = 0.2 * newMean.+ 0.8 * cumMean
      cumMeanRepaired = bounceBackBoundary2(cumMean)

      fn_cum = fn_l(cumMeanRepaired)

      if fn_cum[1] < best_fit[1]
        best_fit = fn_cum
        best_par = cumMeanRepaired
      end

      if fitness[1] <= stopfitness
        msg = "Stop fitness reached."
        break
      end
    end
  end

    log = []

    if log_Ft
       log_Ft = Ft_log
     end
    if log_value
       log_value = value_log[1:iter, :]
     end
    if log_mean
       log_mean = mean_log[1:iter]
     end
    if log_meanCord
       log_meanCord = meanCords_log
     end
    if log_pop
       log_pop = pop_log[:, :, 1:iter]
     end
    if log_bestVal
       log_bestVal = bestVal_log
     end
    if log_worstVal
       log_worstVal = worstVal_log
     end
    if log_eigen
      log_eigen = eigen_log
    end

     res = [
       best_par,
       best_fit,
       convert(Int32, counteval),
       restart_number,
       ifelse(iter >= maxiter, 1, 0),
       msg,
       log,
       iter
     ]

   end
     return res

  end

function rnormFromMatrix(m)
  for i in 1:size(m, 1)
    for j in 1:size(m, 2)
      m[i,j]=randn(m[i, j])
    end
  end
end

function crossprod(x)
    return sum(x.^2)
end

function rastrigin(x)
  d = length(x)
  sum = 10*d
  for i in 1:length(x)
    sum = sum + x[i]*x[i] - 10 * cos(2*pi*x[i])
  end
  return sum
end

function griewank(x)
  sum = 0
  p = 1
  for i in 1:length(x)
    p = p * cos(x[i]/sqrt(i))
    sum = sum + x[i]*x[i]
  end
  sum = sum/4000
  sum += 1
  sum = sum - p
  return sum
end


text:
# default parameters
const default_v_center = 0.0    # m/s
const default_v_step = 250.0    # m/s
const default_v_max = 15.0e3    # m/s
const default_v_width = 410.0   # m/s
const default_v_range_no_mask_change = default_v_max # m/s

include("plan.jl")
export AbstractCCFPlan, BasicCCFPlan
export calc_ccf_v_grid, calc_length_ccf_v_grid

include("calc_ccf.jl")
export ccf_1D, ccf_1D!

include("ccf_template.jl")
export calc_normalized_ccfs, calc_ccf_template

#=
include("convenience/convenience.jl")
export calc_ccf_chunk, calc_ccf_chunk!
export calc_ccf_and_var_chunk, calc_ccf_and_var_chunk!
export calc_ccf_chunklist #, calc_ccf_chunklist!
export calc_ccf_and_var_chunklist #, calc_ccf_and_var_chunklist!
export calc_order_ccfs_chunklist, calc_order_ccfs_chunklist!
export calc_order_ccf_and_vars_chunklist, calc_order_ccf_and_vars_chunklist!
export calc_ccf_chunklist_timeseries
export calc_ccf_and_var_chunklist_timeseries
=#

# Will soon comment out to reduce precompilation time
#include("calc_ccf_old.jl")
#include("convenience_old.jl")
#export ccf_1D_old, ccf_1D_old!

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
"""
# module Find

# Examples

```jldoctest
julia>
```
"""
module Find

using BisectPy: bisect_left, bisect_right

function index(a, x)
    i = bisect_left(a, x)
    (i != length(a) && a[i] == x) ? i : error("Index not found!")
end

end


text:
module TrustPilotCalc

export getTrustScore, playBack

using Dates
using JuliaDB

# references
# https://support.trustpilot.com/hc/en-us/articles/201748946
# http://www.ebc.cat/2015/01/05/how-to-rank-restaurants/
# https://www.evanmiller.org/bayesian-average-ratings.html
# https://medium.com/district-data-labs/computing-a-bayesian-estimate-of-star-rating-means-651496a890ab

"""
TrustPilot review structure
"""
struct TPReview
    author_id::AbstractString
    stars::Float64
    created_at::DateTime
end

"""
TrustPilot score structure
"""
struct TPScore
    score::Float64
    stars::Float64
    date::DateTime
end

"""
  getReviewWeightStatic(review_date, last_review_date)

static variant of review weighting function assuming no decay.
"""
function getReviewWeightStatic(reviewed_at::DateTime, last_reviewed_at::DateTime)
    milli_month = 2592000000
    delta = round((last_reviewed_at - reviewed_at).value / milli_month)
    return if delta < 1 1
        elseif delta < 3 .8
        elseif delta < 6 .5
        elseif delta < 12 .3
        else .1
        end
end

"""
  show(result)

print values of result
"""
function show(result::TPScore)
    printstyled(result.score, " | ", result.stars, " | ", result.date, "\n", color = :blue)
end

"""
  getReviewWeightHL(review_date, last_review_date)

exponential decay (half-life) variant of the `getReviewWeight` function
"""
function getReviewWeightHL(reviewed_at::DateTime, last_reviewed_at::DateTime)
    milli_month = 2592000000 # a month in millis
    initial_weight = 1 # starting weight
    hl = 0.5 # half life
    hlt = 6 # half life t
    delta = round((last_reviewed_at - reviewed_at).value / milli_month)
    weight = initial_weight * exp((log(hl) / hlt) * delta)
    return weight >= 1 ? 1 : weight
end

"""
  getTrustScore(reviews_array)

calculates TrustPilot score given an array of reviews.

returns a `TPScore` object
"""
function getTrustScore(reviews::Array{TPReview})
    first_reviewed_at = mapreduce(x -> x.created_at, min, reviews)
    last_reviewed_at = mapreduce(x -> x.created_at, max, reviews)
    base_reviews = fill(TPReview("seed", 3.5, first_reviewed_at), 7) # avg seed
    all_reviews = [base_reviews; reviews]
    all_weights = mapreduce(x -> getReviewWeightHL(x.created_at, last_reviewed_at), +, all_reviews)
    R = mapreduce(x -> (x.stars * getReviewWeightHL(x.created_at, last_reviewed_at)), +, all_reviews) / all_weights # average user rating for this product
    C = 3.5 # average user rating for all products (seed) aka 3.5
    v = length(reviews)
    m = length(base_reviews)
    w = 1 #v/(v + m) # weight  for R... v/(v + m) v = n rating for this product, m = ... all products aka 7

    S = w * R + (1 - w) * C
    return TPScore(round(S, digits=1), scoreToStars(S), last_reviewed_at)
end

"""
  scoreToStars(trustPilot_score)

rounds TrustPilot score to its equivalent stars value.
"""
function scoreToStars(score::Float64)
    score = round(score, digits=1)
    rem = round(score % .5, digits = 1)
    stars = rem <= .2 ? score - rem : score + (.5 - rem)
    return stars < 1 ? 1 : stars
end

"""
  playBack(csv_file)

imports a TrustPilot reviews CSV export and calculates a TrustPilot score.
"""
function playBack(filename::String; max_yearmon = round(Dates.now(), Dates.Month), schema = (:review_author_id, :review_created_at_utc, :review_star))
    raw_data = loadtable(filename; type_detect_rows = 500)
    working_data = select(raw_data, schema)
    working_data = filter(x -> x.review_created_at_utc < max_yearmon, working_data)
    sort!(working_data, :review_created_at_utc)
    all_reviews = map(x -> TPReview(x.review_author_id, x.review_star, x.review_created_at_utc), working_data)
    return getTrustScore(all_reviews)
end

function playBack(dataset::IndexedTable; max_yearmon = round(Dates.now(), Dates.Month), schema = (:review_author_id, :review_created_at_utc, :review_star))
    raw_data = dataset
    working_data = select(raw_data, schema)
    working_data = filter(x -> x.review_created_at_utc < max_yearmon, working_data)
    sort!(working_data, :review_created_at_utc)
    all_reviews = map(x -> TPReview(x.review_author_id, x.review_star, x.review_created_at_utc), working_data)
    return getTrustScore(all_reviews)
end

end # module


text:
# This file includes an example of memory operation.
# In this example `mem` delays its input for one step size time.

using Causal
using Plots

# Construct a model
ti, dt, tf = 0., 1., 100.
numtaps = 5  # Number of buffer taps in Memory.
delay = dt   # One step size delay.
model = Model(clock=Clock(ti, dt, tf))

@defmodel model begin
    @nodes begin
        gen = RampGenerator()
        mem = Memory(delay=dt, numtaps=numtaps)
        writer = Writer(input=Inport(2))
    end
    @branches begin
        gen     =>  mem
        mem[1]  =>  writer[1]
        gen[1]  =>  writer[2]
    end
end

# Simulate model
simulate!(model, ti, dt, tf)

# Read simulation data
t, x = read(getnode(model, :writer).component)
u = getnode(model, :gen).component.readout.(t .- delay)
err = u - x[:, 2]

# Plots simulation data
n1, n2 = 1, 5
marker = (:circle, 2)
p1 = plot(t[n1:n2], x[n1:n2, 1], label=:mem, marker=marker)
plot!(t[n1:n2], x[n1:n2, 2], label=:gen, marker=marker)
plot!(t[n1:n2], u[n1:n2], label=:real, marker=marker)
p2 = plot(t[n1:n2], err[n1:n2], label=:error, marker=marker)
display(plot(p1, p2, layout=(2,1)))


text:
function random_neighborhood(mapf::MAPF, neighborhood_size::Integer)
    return sample(1:nb_agents(mapf), neighborhood_size; replace=false)
end

function random_neighborhood_collision_degree(
    solution::Solution, mapf::MAPF, neighborhood_size::Integer
)
    wv = StatsBase.weights(collision_degrees(solution, mapf) .+ 1)
    return sample(1:nb_agents(mapf), wv, neighborhood_size; replace=false)
end

function conflicting_neighborhood(solution::Solution, mapf::MAPF)
    A = nb_agents(mapf)
    neighborhood = Set{Int}()
    for a in 1:A
        for b in 1:(a - 1)
            if have_conflict(a, b, solution, mapf)
                push!(neighborhood, a)
                push!(neighborhood, b)
                break
            end
        end
    end
    return collect(neighborhood)
end

function remove_agents!(solution::Solution, agents)
    backup = Dict(a => solution[a] for a in agents)
    for a in agents
        solution[a] = Path()
    end
    return backup
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(DWARF.LineTableSupport.decode_next, (Base.AbstractIOBuffer{Array{UInt8, 1}}, DWARF.LineTableSupport.Header{UInt32},))
    precompile(DWARF.Attributes.readorskip, (DWARF.LightDIERef, UInt8, Symbol, Base.Val{:skip},))
    precompile(DWARF.Attributes.readorskip, (DWARF.LightDIERef, UInt8, Symbol, Base.Val{:read},))
    precompile(DWARF.LineTableSupport.step, (Base.AbstractIOBuffer{Array{UInt8, 1}}, DWARF.LineTableSupport.StateMachine,))
    precompile(DWARF.skip_attributes, (DWARF.LightDIERef, DWARF.AbbrevTableEntry,))
    precompile(DWARF.CallFrameInfo.evaluage_op, (DWARF.CallFrameInfo.RegStates, Base.AbstractIOBuffer{Array{UInt8, 1}}, DWARF.CallFrameInfo.CIE,))
    precompile(DWARF.LineTableSupport.pc_adv!, (DWARF.LineTableSupport.StateMachine, Int64,))
    precompile(DWARF.CallFrameInfo.realize, (DWARF.CallFrameInfo.CIERef,))
    precompile(DWARF._searchforip, (DWARF.ChildIterator, UInt64,))
    precompile(DWARF.LineTableSupport.Type, (Type{DWARF.LineTableSupport.RegisterState}, Bool,))
    precompile(DWARF.decode, (Array{UInt8, 1}, Int64, Type{DWARF.SLEB128}, DataType,))
    precompile(DWARF.LineTableSupport.pc_adv!, (DWARF.LineTableSupport.StateMachine, Base.GMP.BigInt,))
    precompile(DWARF.next, (DWARF.UnitIterator{DWARF.DWARFCUHeader}, UInt32,))
    precompile(DWARF.next, (DWARF.UnitIterator{DWARF.DWARFCUHeader}, Int64,))
    precompile(DWARF.LineTableSupport.read, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.LineTableSupport.FileEntry},))
    precompile(DWARF.read, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.SLEB128},))
    precompile(DWARF.LineTableSupport.readstring, (Base.AbstractIOBuffer{Array{UInt8, 1}},))
    precompile(DWARF.extract_attribute, (DWARF.LightDIERef, UInt8,))
    precompile(DWARF.Expressions.evaluate_generic, (DWARF.Expressions.StateMachine{UInt64}, Array{UInt8, 1}, Any, Any, Any, Symbol,))
    precompile(DWARF.size_to_inttype, (UInt8,))
    precompile(DWARF.read, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.DWARFCUHeader}, Symbol,))
    precompile(DWARF.skip_children, (Base.AbstractIOBuffer{Array{UInt8, 1}}, DWARF.LightDIERef,))
    precompile(DWARF.extract_attribute, (DWARF.LightDIERef, UInt8, Symbol,))
    precompile(DWARF.LineTableSupport.unpack, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.LineTableSupport.HeaderStub},))
    precompile(DWARF.Attributes.readorskip, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.ULEB128}, Symbol, Base.Val{:skip},))
    precompile(DWARF.pc_range, (DWARF.LightDIERef,))
    precompile(DWARF.LineTableSupport.Type, (Type{DWARF.LineTableSupport.HeaderStub{UInt32}}, UInt32, UInt16, UInt32, UInt8, Int64, UInt8, Int8, Int8, UInt8,))
    precompile(DWARF._searchforip, (DWARF.UnitIterator{DWARF.DWARFCUHeader}, UInt64,))
    precompile(DWARF.LineTableSupport.pc_adv!, (DWARF.LineTableSupport.StateMachine, UInt8,))
    precompile(DWARF.read, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.AbbrevTableEntry}, Symbol,))
    precompile(DWARF.Expressions.evaluate_simple_location, (DWARF.Expressions.StateMachine{UInt64}, Array{UInt8, 1}, Any, Any, Any, Symbol,))
    precompile(DWARF.LineTableSupport.read_header, (Base.AbstractIOBuffer{Array{UInt8, 1}},))
    precompile(DWARF.LineTableSupport.start, (DWARF.LineTableSupport.LineTable,))
    precompile(DWARF.Expressions.Type, (Type{DWARF.Expressions.StateMachine{UInt64}},))
    precompile(DWARF.LineTableSupport.Type, (Type{DWARF.LineTableSupport.Header{UInt32}}, DWARF.LineTableSupport.HeaderStub{UInt32}, Array{UInt8, 1}, Array{String, 1}, Array{DWARF.LineTableSupport.FileEntry, 1},))
    precompile(DWARF.LineTableSupport.header_type, (DWARF.LineTableSupport.Header{UInt32},))
    precompile(DWARF.CallFrameInfo.encoding_type, (UInt8,))
    precompile(DWARF.Attributes.convert, (Type{UInt64}, DWARF.Attributes.SectionReference,))
    precompile(DWARF.Attributes.Type, (Type{DWARF.Attributes.SectionReference}, UInt32,))
    precompile(DWARF.Type, (Type{DWARF.AttributeSpecification}, UInt64, UInt8,))
    precompile(DWARF.CallFrameInfo.getindex, (DWARF.CallFrameInfo.RegStates, UInt8,))
    precompile(DWARF.Attributes.Type, (Type{DWARF.Attributes.StrTableReference}, UInt32,))
    precompile(DWARF.Type, (Type{DWARF.LightDIERef}, Base.AbstractIOBuffer{Array{UInt8, 1}}, DWARF.CU, Int64,))
    precompile(DWARF.Attributes.Type, (Type{DWARF.Attributes.DwarfExpr{UInt64}}, Void,))
    precompile(DWARF.Attributes.Type, (Type{DWARF.Attributes.DwarfExpr{UInt64}}, Array{UInt8, 1},))
    precompile(DWARF.start, (DWARF.ChildIterator,))
    precompile(DWARF.seekentry, (DWARF.LightDIERef,))
    precompile(DWARF.read, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.ULEB128},))
    precompile(DWARF.LineTableSupport.header_type, (DWARF.LineTableSupport.HeaderStub{UInt32},))
    precompile(DWARF.CallFrameInfo.convert, (Type{Int64}, DWARF.CallFrameInfo.RegNum,))
    precompile(DWARF.CallFrameInfo.Type, (Type{DWARF.CallFrameInfo.RegStates},))
    precompile(DWARF.CallFrameInfo.setindex!, (DWARF.CallFrameInfo.RegStates, DWARF.CallFrameInfo.Offset, DWARF.CallFrameInfo.RegNum,))
    precompile(DWARF.Type, (Type{DWARF.ULEB128}, Base.GMP.BigInt,))
    precompile(DWARF.CallFrameInfo.fde_range, (DWARF.CallFrameInfo.FDERef, DWARF.CallFrameInfo.CIE,))
    precompile(DWARF.LineTableSupport.state_step, (Base.AbstractIOBuffer{Array{UInt8, 1}}, DWARF.LineTableSupport.StateMachine,))
    precompile(DWARF.CallFrameInfo.evaluate_program, (Base.AbstractIOBuffer{Array{UInt8, 1}}, UInt64, DWARF.CallFrameInfo.CIE, DWARF.CallFrameInfo.RegStates, UInt64,))
    precompile(DWARF.read, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Type{DWARF.AbbrevTableSet}, Symbol,))
    precompile(DWARF.readnthattr, (DWARF.LightDIERef, DWARF.AbbrevTableEntry, Symbol, Int64,))
    precompile(DWARF.convert, (Type{UInt64}, DWARF.ULEB128,))
end


text:
"""

    const TileMap{O} = ObjectIndexableArray{Bool, 3, BitArray{3}, O}

An instance of `TileMap`, referred to as `tile_map` here, wraps an `array` of type `BitArray{3}` and is of size `(num_objects, height, width)`, which encodes information about the presence or absence of objects across the tiles using Boolean values. Each tile can contain multiple objects, which is captured by a multi-hot encoding along the first dimension (`num_objects` dimension) of the `array`.

"""
const TileMap{O} = ObjectIndexableArray{Bool, 3, BitArray{3}, O}

TileMap(objects::Tuple{Vararg{AbstractObject}}, grid::BitArray{3}) = TileMap{typeof(objects)}(grid)

TileMap(objects::Tuple{Vararg{AbstractObject}}, height::Integer, width::Integer) = TileMap(objects, falses(length(objects), height, width))

get_num_objects(tile_map::TileMap) = size(tile_map, 1)
get_height(tile_map::TileMap) = size(tile_map, 2)
get_width(tile_map::TileMap) = size(tile_map, 3)


text:
using FluidInducedAseismicSlip
using Plots

function test_slip()
    (x, δ, λ) = injection_analytical_gs(0.5, 500)

    # Analytical slip
    x_a = collect(0:0.05:1.0)
    δ_a = [0.23366202, 0.23081724, 0.22426165, 0.21513089, 0.20404685, 0.19145812, 0.17772295, 0.16314401, 0.14798653, 0.13248893, 0.11686999, 0.101334019, 0.086075209, 0.071281873, 0.057141201, 0.043845559, 0.031602114, 0.020649782, 0.0112944788, 0.0040083095, 0.0]

    # Plot
    plot(x, δ, label="This package")
    plot!(x_a, δ_a, label="Viesca (2021)")
end

test_slip()

max_stars_repo_name: ajacquey/FluidInducedAseismicSlip

max_issues_repo_name: ajacquey/FluidInducedAseismicSlip

max_forks_repo_name: ajacquey/FluidInducedAseismicSlip


text:
using SparseArrays

include("../hypergraphs.jl")
include("../SparseCard.jl")
include("../pwl_approx.jl")

## The HyperModualrity package includes a number of hypergraph datasets we can use for testing
using HyperModularity
hypermodularity_datasets()

## Load dataset and solve
dataset = "contact-primary-school-classes"
# dataset = "mathoverflow-answers"
maxsize = 25	# max hyperedge size
minsize = 2	# min hyperedge size
return_labels = true
Edges, Labels = read_hypergraph_data(dataset,maxsize,minsize,return_labels)
pb = 0.2:0.1:0.3
EdgesW = RandomSplittingFunctions(Edges,pb)
P = sortperm(Labels)
n = length(Labels)
refined = false
epsilon = 0.1*ones(maxsize)

## Solve the problem with two different max-flow subroutines
eS_pr, cutval_pr, approx_pr, flowtime_pr, reduce_time, worst_eps = SparseCard(Edges,EdgesW,n,epsilon;refined = false)
eS_gur, cutval_gur, approx_gur,flowtime_gur, reduce_time, worst_eps = SparseCard(Edges,EdgesW,n,epsilon;flowsolver = "gurobi")

## print output details

println("Reduce time = $reduce_time \t Approx =  $worst_eps")
println("PushRelabel: $cutval_pr \t  $approx_pr \t $flowtime_pr \t")
println("Gurobi-Flow: $cutval_pr \t  $approx_pr \t $flowtime_gur")

## More detail for reduction
A, svec, tvec = CardDSFM_reduction(Edges,EdgesW,n,epsilon,refined)
N = size(A,1)

## Solve the max-flow with push-relabel code
F = maxflow(A,svec,tvec,0)
Sall = F.source_nodes[2:end].-1
ct = F.cutvalue
Small = source_nodes_min(F,1e-10)[2:end].-1

S = intersect(Sall,1:n)
eS = zeros(n)
eS[S] .= 1
cutS = eval_fS(Edges,EdgesW,eS)

## Find the best cut with certain nodes fixed
T = setdiff(1:n,S)
S2, cutval = eval_min_dircut(A,svec,tvec,S,T)

## Solve the max-flow with Gurobi
Svec, flowval, f = Gubori_MaxFlow(A,svec,tvec,1)
eS = Svec[1:n]
cut_gurobi = eval_fS(Edges,EdgesW,Svec)

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
import ..GroupSlices:
    groupslices,
    groupinds

import ..Embeddings:
    AbstractEmbedding

"""
    marginal_visitation_freq(along_which_axes::Union{Int, Vector{Int}, AbstractUnitRange{Int}},
        visited_bins::Array{T, 2}) where T

Calculate marginal or joint visitation frequencies for a pre-binned set of
points. Each column in `visited_bins` corresponds to a unique point, and
contains the coordinate or integer representation of the bin containing that
point.

## Arguments
- `visited_bins`: Contains a representation of the boxes containing each point.
- `along_which_axes`: Controls which axes to take the marginal visitation
    frequencies along. Setting `along_which_axes` to a range 1:D, where
    D is the dimension of the corresponding state space, corresponds to taking the
    joint visitation frequency (basically, a D-dimensional histogram).
- `npts`: The number of points in the embedding, which is used as a
    normalization factor.
"""
function marginal_visitation_freq(
        along_which_axes::Union{Int, Vector{Int}, AbstractUnitRange{Int}},
        visited_bins::AbstractArray{T, 2}) where T
    if typeof(along_which_axes) <: Int
        along_which_axes = [along_which_axes]
    end
    # The number of points to which we have assigned bins
    npts = size(visited_bins, 2)

    # Assuming there are M unique bins (or, equivalently, integer tuples),
    # groupslices(visited_bins_inds, 2) returns a vector of indices,
    # in which each visited bin label is associated with the i-th unique
    # bin. For example, if slices = [1, 2, 3, 1, 2, 3, 7, ...], then points
    # p₁ and p₄ are repeated (have identical labels), p₂ and p₅ are repeated,
    # and p₃ and p₆ are repeated.
    #
    # Subsetting some of the rows are equivalent to estimating the marginal
    # along the corresponding axes. If along_which_axes = 1:n_axes, then we're
    # computing the joint frequency distribution.
    slices = groupslices(visited_bins[along_which_axes, :], 2)

    # The groupinds function takes these indices and groups the indices
    # corresponding to repeated points. For the example in the comment
    # above, group_repeated_inds = [[1, 4], [2, 5], [3, 6],
    # [7, indices of potential points shared with p₇], ....].
    group_repeated_inds = groupinds(slices)

    # Computing the joint probability for pᵢ is now just a matter of counting
    # the number of times each visited bin is visited. The number of times
    # the ith visited bin is visited by the orbit is then simply counting
    # how many elements there are in group_repeated_inds[i].
    n_visited_states = length(group_repeated_inds)

    # We'll now loop over one bit at a time, counting how many times
    # the bin is visited, then get visitation frequency for that bin
    # by dividing the number of times it is visited by the total
    # number of points.
    m = Vector{Float64}(undef, n_visited_states) # preallocate marginal vector
    @inbounds for i = 1:n_visited_states
        m[i] = length(group_repeated_inds[i]) / npts
    end

    return m
end

"""
    marginal_visitation_freq(
        along_which_axes::Union{Int, Vector{Int}, AbstractUnitRange{Int}},
        points::Array{T, 2},
        ϵ) where T

Compute marginal visitation frequencies for a set of `points`, given a partition
scheme `ϵ`. The following `ϵ` will work:

* `ϵ::Int` divide each axis into `ϵ` intervals of the same size.
* `ϵ::Float` divide each axis into intervals of size `ϵ`.
* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size.
* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`.

The argument `along_which_axes` controls which axes to take the marginal
visitation frequencies along. Setting `along_which_axes` to a range 1:D,
where D is the dimension of the corresponding state space, corresponds to
taking the joint visitation frequency.
"""
function marginal_visitation_freq(
            along_which_axes::Union{Int, Vector{Int}, AbstractUnitRange{Int}},
            points::AbstractArray{T, 2},
            ϵ) where T
    # Make sure that the array contains points as columns.
    if size(points, 1) > size(points, 2)
        error("The dimension of the dataset exceeds the number of points.")
    end

    N = size(points, 1)
    visited_bins = assign_bin_labels(points, ϵ)
    marginal_visitation_freq(along_which_axes, visited_bins)
end

"""
    marginal_visitation_freq(
        along_which_axes::Union{Int, Vector{Int}, AbstractUnitRange{Int}},
        E::AbstractEmbedding,
        ϵ)

Compute marginal visitation frequencies for an embedding, given a partition
scheme `ϵ`. The following `ϵ` will work:

* `ϵ::Int` divide each axis into `ϵ` intervals of the same size.
* `ϵ::Float` divide each axis into intervals of size `ϵ`.
* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size.
* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`.

The argument `along_which_axes` controls which axes to take the marginal
visitation frequencies along. Setting `along_which_axes` to a range 1:D,
where D is the dimension of the corresponding state space, corresponds to
taking the joint visitation frequency.
"""
function marginal_visitation_freq(
            along_which_axes::Union{Int, Vector{Int}, AbstractUnitRange{Int}},
            E::AbstractEmbedding, ϵ)
    marginal_visitation_freq(along_which_axes, E.points, ϵ)
end

marginal_visitfreq = marginal_visitation_freq

export
marginal_visitation_freq


text:
# ---
# title: 731. My Calendar II
# id: problem731
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Ordered Map
# link: <https://leetcode.com/problems/my-calendar-ii/description/>
# hidden: true
# ---
#
# Implement a `MyCalendarTwo` class to store your events. A new event can be
# added if adding the event will not cause a **triple** booking.
#
# Your class will have one method, `book(int start, int end)`. Formally, this
# represents a booking on the half open interval `[start, end)`, the range of
# real numbers `x` such that `start <= x < end`.
#
# A _triple booking_ happens when **three** events have some non-empty
# intersection (ie., there is some time that is common to all 3 events.)
#
# For each call to the method `MyCalendar.book`, return `true` if the event can
# be added to the calendar successfully without causing a **triple** booking.
# Otherwise, return `false` and do not add the event to the calendar.
#
# Your class will be called like this: `MyCalendar cal = new MyCalendar();`
# `MyCalendar.book(start, end)`
#
# **Example 1:**
#
#
#
#     MyCalendar();
#     MyCalendar.book(10, 20); // returns true
#     MyCalendar.book(50, 60); // returns true
#     MyCalendar.book(10, 40); // returns true
#     MyCalendar.book(5, 15); // returns false
#     MyCalendar.book(5, 10); // returns true
#     MyCalendar.book(25, 55); // returns true
#     Explanation:
#     The first two events can be booked.  The third event can be double booked.
#     The fourth event (5, 15) can't be booked, because it would result in a triple booking.
#     The fifth event (5, 10) can be booked, as it does not use time 10 which is already double booked.
#     The sixth event (25, 55) can be booked, as the time in [25, 40) will be double booked with the third event;
#     the time [40, 50) will be single booked, and the time [50, 55) will be double booked with the second event.
#
#
#
#
# **Note:**
#
#   * The number of calls to `MyCalendar.book` per test case will be at most `1000`.
#   * In calls to `MyCalendar.book(start, end)`, `start` and `end` are integers in the range `[0, 10^9]`.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module CVModule # All functions/types needed for closed-loop CV system solver

importall NumericalIntegration
importall MAT
importall Interpolations
importall CSV
importall ZMQ
importall Missings

# conversion factors and solver parameters
include("conversions.jl")
include("solverparams.jl")

# function for loading 1D domain empirical data from CSV file
include("loadtexttree.jl")

# type definitions
include("buildbranches.jl")
include("buildheart.jl")
include("buildvenacava.jl")
include("buildlungs.jl")
include("buildliver.jl")
include("buildcns.jl")
include("buildhemo.jl")
include("buildall.jl")

# memory allocators for solution variables
include("calcbranchprops.jl")
include("discretizebranches.jl")
include("assignterminals.jl")
include("discretizeperiphery.jl")
include("discretizeheart.jl")
include("discretizelungs.jl")
include("discretizeliver.jl")
include("discretizecns.jl")
include("updatediscretization.jl")

# initial conditions for each submodel
include("applybranchics.jl")
include("applyperipheryics.jl")
include("applyheartics.jl")
include("applylungics.jl")
include("applyliverics.jl")
include("applycnsics.jl")
include("applycustomics.jl")

# heart model
include("elastancefn.jl")
include("elastancemodel.jl")
include("setnumbeats.jl")

# volume tracker, error corrector for mass conservation
include("updatevolumes.jl")
include("correctvolume.jl")

# Lax-Wendroff method for interior points of arterial network
include("predictorfluxes.jl")
include("predictorstep.jl")
include("correctorfluxes.jl")
include("correctorstep.jl")

# 0D-1D coupling and 0D updates
include("applyendbcs.jl")
include("coupledistal.jl")
include("endinvariants.jl")
include("updateterms.jl")
include("updateliver.jl")
include("updatevc.jl")
include("updaterh.jl")
include("updatelungs.jl")
include("updatela.jl")
include("coupleproximal.jl")
include("rootinvariant.jl")
include("fav.jl")
include("Jav.jl")
include("newtonav.jl")
include("linesearch.jl")
include("fdist.jl")
include("Jdist.jl")
include("newtondist.jl")
include("linedist.jl")

# 1D interior junction updates
include("splitinvariants.jl")
include("solvesplits.jl")
include("newton.jl")
include("fsingle.jl")
include("fdouble.jl")
include("ftriple.jl")
include("fquad.jl")
include("Jsingle.jl")
include("Jdouble.jl")
include("Jtriple.jl")
include("Jquad.jl")

# 1D pressure updates
include("arterialpressure.jl")

# regulation
include("regulateall.jl")
include("regulateperiphery.jl")
include("reflexpressure.jl")
include("cnsactivations.jl")
include("regulateheart.jl")

# systemic arterial hemorrhage/tourniquet application
include("applytourniquet.jl")
include("applyhemoics.jl")
include("model1dhemo.jl")

# ZMQ message passing
include("senddata.jl")
include("sc.jl")

export CVSystem
export Heart
export ArterialBranches
export SolverParams
export VenaCava
export Lungs
export Liver
export CNS
export Hemorrhage

export loadtexttree
export calcbranchprops!
export discretizebranches!
export assignterminals!
export discretizeperiphery!
export discretizeheart!
export discretizelungs!
export discretizeliver!
export discretizecns!
export applybranchics!
export applyperipheryics!
export applyheartics!
export applylungics!
export applyliverics!
export applycnsics!
export applycustomics!
export elastancefn!
export elastancemodel!
export updatevolumes!
export buildall
export predictorfluxes!
export predictorstep!
export correctorfluxes!
export correctorstep!
export applyendbcs!
export coupledistal!
export endinvariants!
export newtondist!
export linedist!
export fdist
export Jdist
export updateterms!
export updateliver!
export updatevc!
export updaterh!
export updatelungs!
export updatela!
export coupleproximal!
export rootinvariant!
export newtonav!
export fav
export Jav
export linesearch
export splitinvariants!
export solvesplits!
export newton!
export fsingle
export fdouble
export ftriple
export fquad
export Jsingle
export Jdouble
export Jtriple
export Jquad
export arterialpressure!
export regulateall!
export regulateperiphery!
export setnumbeats!
export reflexpressure!
export cnsactivations!
export regulateheart!
export updatediscretization!
export applytourniquet!
export modelhemo!
export applyhemoics!
export senddata
export sc

end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
#Setup

include("make_plots_distributed.jl")

using .MakePlots

println("Ready to set up workers.")
flush(stdout)
MakePlots.setup() # setup workers based on cpu counts

#Run simulation

d_tasks = parse(Int, ENV["SLURM_ARRAY_TASK_ID"]) # 1 - 3

job_id = "rz_0.001"
samples = 1000*20000 # 20,000 samples take about 1min per cpu (1 plot, 1 dist, 10 angles)
dists = [(1 + 2 * d_tasks):2:(2 + 2 * d_tasks)...] # code distance: 3 or 5 or 7 separately
plots = [4] # list of plot types
plot_id = 4 # type of plots to be made

println("Simulation begins... (", samples, " samples)")
flush(stdout)

results = MakePlots.dist_calc_all(job_id, dists, samples, plots) # start computing and wait until complete
println("Simulation completed.")

collected = MakePlots.fetch_finished(results, job_id, plot_id)
println(collected)

max_stars_repo_name: yongshanding/VLQ

max_issues_repo_name: yongshanding/VLQ

max_forks_repo_name: yongshanding/VLQ


text:
export create_∂

# Assumption: we don't calculate derivatives for interpolated fields.  In other words, we
# calculate derivatives for the original fields defined on Yee's grid.  (There are
# exceptions; see below.)  This assumption simplifies implementation of the symmetry
# boundary condition in the derivative operators greatly.
#
# For example, suppose we use the primal grid planes as the E-field planes.  Then, for the
# symmetry boundary condition (which is the PEC boundary condition in this case), the fields
# defined on the domain boundaries are the tangential  E-fields and normal H-fields.  Both
# fields are zero on the symmetry boundaries.  Therefore, if we assume that we calculate
# derivatives only for those fields defined on Yee's grid, we do not have to construct
# separate derivative operators for the tangential E-field and normal H-field.  (Note that
# the derivative operator for the tangential E-field is used in curl for E, and the
# derivative operator for the normal H-field is used in divergence for H.)
#
# If we allow application of derivative operators to interpolated fields, there are cases
# where the boundary condition for those fields cannot be implemented in the derivative
# operators simply by zeroing the interpolated fileds.  For instance, in the aformentioned
# system, suppose we interpolate the normal E-field on the boundary.  Then, because the
# normal E-fields remain the same around the boundary, the interpolation on the boundary
# surface does not lead to zero.  Therefore, in the derivative operator for the normal
# E-field, we cannot implement the symmetry boundary condition simply by zeroing the fields.
# Supporting such interpolated fields would make the code for the derivative operators very
# convoluted, so we will assume that we don't deal with such cases.
#
# There are exceptions for this assumption.  As long as the interpolated fields on the
# symmetry boundaries are zero, our derivative operator code can be used to differentiate
# such interpolated fields.  For example, again in the aformentioned system, even though the
# normal D-field is not zero on the symmetry boundary, a portion of it could be zero in the
# following situation.  Suppose we are dealing with anisotropic materials.  Then, if we are
# dealing with the z-normal boundaries, Dz is the normal D-field but it is expressed as
# Dz = εzx Ex + εzy Ey + εzz Ez = Dzx + Dzy + Dzz.  Here, even though Dzz is nonzero because
# Ez is nonzero, Dzx and Dzy are zero because Ex and Ey are tangential to the PEC boundary.
# Therefore, we can use our derivative operator to differentiate Dzx and Dzy along the
# z-direction.  (I don't think there will be cases where we would need to differentiate Dzx
# and Dzy, but I mention this for consistency with the averaging operators in mean.jl,
# because we do need to take averages for Dzx and Dzy.)

## Difference operators ##

# Wrapper to create the discrete partial derivative by default
create_∂(nw::Integer,  # 1|2|3 for x|y|z; 1|2 for horizontal|vertical
         isfwd::Bool,  # true|false for forward|backward difference
         N::AbsVecInteger,  # size of grid
         ∆w⁻¹::Number=1.0,  # inverse of spatial discretization
         isbloch::Bool=true,  # boundary condition in w-direction
         e⁻ⁱᵏᴸ::Number=1.0) =  # Bloch phase factor
    create_∂(nw, isfwd, N, fill(∆w⁻¹, N[nw]), isbloch, e⁻ⁱᵏᴸ)  # fill: create vector of ∆w⁻¹

# Wrapper to convert AbstractVector's to SVec's
create_∂(nw::Integer,  # 1|2|3 for x|y|z; 1|2 for horizontal|vertical
         isfwd::Bool,  # true|false for forward|backward difference
         N::AbsVecInteger,  # size of grid
         ∆w⁻¹::AbsVecNumber,  # inverse of spatial discretization; vector of length N[nw]
         isbloch::Bool=true,  # boundary condition in w-direction
         e⁻ⁱᵏᴸ::Number=1.0) =  # Bloch phase factor
    (K = length(N); M = prod(N); dropzeros!(sparse(create_∂info(nw, isfwd, SInt{K}(N), ∆w⁻¹, isbloch, e⁻ⁱᵏᴸ)..., M, M)))

# I need to figure out whether the ±1 entries of the backward difference operator is always
# the transpose of the forward difference operator for all boundary conditions. (∆w division
# factors vary, though.)  This was the case in the MATLAB code, but in the Julia code I
# changed the treatment of PDC, so let's make sure about this again.
#
# For the symmetry boundary for the E-field (then the E-fields are on the primal grid planes),
# assuming E that is correctly zeroed at the negative boundary is supplied to the forward
# difference operator, we can take the ±1 pattern of the forward difference operator for the
# Bloch boundary as the ±1 pattern of the forward difference operator for the symmetry
# boundary, because E = 0 is correctly used for the value at the positive boundary because
# the field value is copied from the negative boundary.
#
# Now, let's think about the backward difference operator for H.  Because the boundary is
# the symmetry boundary for E, the ghost H₀, which is before the negative boundary, must be
# the same as the non-ghost H₁.  Therefore, this leads to the first difference being H₁-H₀ =
# H₁-H₁ = 0, which means that the first row of the backward difference operator must be
# empty.  However, when the forward difference operator for the symmetry boundary is created
# the same as that for Bloch, its transpose does not have an empty first row!
#
# For the backward differce operater to be the transpose of the forward difference operator,
# I need to create the forward difference operator such that E₁ is zeroed.  This turns out
# to work.  See my notes on Sep/06/2017 in RN - MaxwellFDM.jl.

# Creates the w-directional difference matrix, with division by ∆w's.
function create_∂info(nw::Integer,  # 1|2|3 for x|y|z; 1|2 for horizontal|vertical
                      isfwd::Bool,  # true|false for forward|backward difference
                      N::SVec{K,Int},  # size of grid
                      ∆w⁻¹::AbsVecNumber,  # spatial discretization; vector of length N[nw]
                      isbloch::Bool,  # boundary condition in w-direction
                      e⁻ⁱᵏᴸ::Number  # Bloch phase factor
                      ) where {K}
    M = prod(N)
    Nw = N[nw]
    ŵ = SVec(ntuple(identity,Val(K))) .== nw  # unit vector in w-direction; [0,true,0] for w == y
    ns = isfwd ? 1.0 : -1.0  # number for sign

    # Below, when constructing I, J's, V's, note that a tuple of array subscripts (i,j,k)
    # indicates a row index of the final matrix.  In other words, the entries with the same
    # subscripts correspond to the same *output* field.  So, when setting I[i,j,k], J₀[i,j,k],
    # Jₛ[i,j,k], V₀[i,j,k], Vₛ[i,j,k], we must know that they are for the entries in the
    # same row index (or the same output field).  Specifically,
    #
    # - I is the row index itself, so it is the identity map: I[i,j,k] = LinearIndices(N)[i,j,k].
    #
    # - J₀[i,j,k] is the column index of the diagonal entry in the row subscripted (i,j,k).
    # Because the row and column indices are the same for the diagonal entries, J₀ = I and
    # not created separately.
    #
    # - Jₛ[i,j,k] is the column index of the off-diagonal entry in the row subscripted
    # (i,j,k).
    #
    # - V₀[i,j,k] is the value of the diagonal entry in the row subscripted (i,j,k).
    #
    # - Vₛ[i,j,k] is the value fo the off-diagonal entry in the row subscripted (i,j,k).

    I = reshape(collect(1:M), N.data)

    Jₛ = reshape(collect(1:M), N.data)
    shifts = -ns * ŵ  # shift vector for sign; [0,-1,0] for w == y and ns = +1
    Jₛ = circshift(Jₛ, shifts.data)

    # Align ∆w⁻¹ in the w-direction.
    vec1 =  @SVector ones(Int,K)
    sizew = @. !ŵ * vec1 + ŵ * N  # [1,Ny,1] for w == y
    ∆W⁻¹ = reshape(∆w⁻¹, sizew.data)

    # Construct the values of the diagonal and off-diagonal nonzero entries of the matrix
    # for the periodic (≠ Bloch) boundary condition.  (These values will be updated for
    # specific boundary conditions later.)
    #
    # To figure out the entries of V₀ and Vₛ, it is most convenient to figure out the entire
    # 3D arrays V₀ and Vₛ, instead of entry by entry.  In doing so, note that the entries
    # with the same subscripts correspond to the same *output* field (= same row index).
    T = promote_type(eltype(∆w⁻¹), eltype(e⁻ⁱᵏᴸ))
    V₀ = -ns .* ones(T, N.data) .* ∆W⁻¹  # values of diagonal entries
    Vₛ = ns .* ones(T, N.data) .* ∆W⁻¹  # values of off-diagonal entries

    # Modify V₀ and Vₛ according to the boundary condition.  What we do is basically to take
    # the operator for the periodic boundary condition (constructed above) as a template and
    # then modify it for the actual boundary condition, either by applying the Bloch phase
    # factors for the Bloch boundary condition or by zeroing some diagonal and off-diagonal
    # entries for the symmetry boundary condition.  Because this modification is due to the
    # boundary condition, the modification needs to be made only at the array entries with
    # the first or last subscript of the Cartesion direction in which the differentiation is
    # taken.
    #
    # The outline of the procedure is as follows:
    #
    # 1. Detertmine the rows in which the modification should be made.  This is done by
    # figuring out the array subscripts (i,j,k) of the output derivative array entries to
    # which the input boundary fields contribute.
    #
    # 2. In each of these rows, determine the column in which the modification should be
    # made.  Because each row has only two nonzero entries (one diagoal and one off-diagonal),
    # this step is to determine whether the modification should be made at the diagonal
    # entry or the off-diagonal entry in the row currently concerned.  If the array
    # subscripts of the input fields contributing to the currently concerned output
    # derivatives are the same as the array subscripts of the output derivatives, the
    # modification is made in V₀ (diagonal).  If they are different, the modification is
    # made in Vₛ (off-diagonal).
    #
    # 3. In the determined rows, modify either V₀ or Vₛ depending on the above step.  This
    # is to modify V₀[i,j,k] or Vₛ[i,j,k], where (i,j,k) is the array subscripts specifying
    # the determined row.
    #
    #
    # For the final sparsity patterns of the operators, see my notes on September 6, 2017 in
    # RN - MaxwellFDM.jl.nb
    #
    # Below, V[Base.setindex(axes(V), iw:iw, nw)...] mimics the implementation of slicedim
    # and means V[:,iw:iw,:] for w = y.  (axes(V) = (1:Nx,1:Ny,1:Nz) and setindex(t,v,i)
    # does t[i] = v, so Base.setindex(axes(V),iw:iw,nw) creates (1:Nx,iw:iw,1:Nz) for w = y.)
    # The use of iw:iw instead of iw is to support 1D arrays.  If V is 1D, then V[iw] is a
    # scalar and the dot equal on V[iw] fails.
    if isbloch
        # Application of the aformentioned procedure:
        #
        # 0. Goal
        # Our goal is to create the ghost fields on one boundary by multiplying the Bloch
        # phase factors to the nonghost fields on the opposite boundary.  Note that w is the
        # direction of differentiation.
        #
        # 1. Determination of the indices of the rows to modify
        # For isfwd = true, we are taking the forward difference between primal grid planes.
        # The ghost input fields are at the positive-end boundary, and they contritube to
        # the derivatives at the last dual grid points, so the rows to modify are
        # subscripted (i,j,k) with iw ≡ (i,j,k)[w] = Nw.
        #
        # For isfwd = false, we are taking the backward difference between dual grid planes.
        # The ghost input fields are before the negative-end boundary, and they contribute
        # to the derivatives at the first primal grid points, so the rows to modify are
        # subscripted (i,j,k) with iw ≡ (i,j,k)[w] = 1.
        #
        # 2. Determination of the column in each row to modify
        # For isfwd = true, the w-subscript of the ghost input field is greater by 1 than
        # that of the output derivative.
        #
        # For isfwd = false, the w-subscript of the ghost input field is less by 1 than that
        # of the output derivative.
        #
        # The nonghost input fields bing brought to these ghost fields (from the opposite
        # boundary) are different from the input fields subscripted the same as the output
        # derivatives.  Therefore, the matrix entries to modify by multiplying the Bloch
        # phase factors with are off-diagonal entries in Vₛ
        #
        # 3. Modification of the values in V
        # For isfwd = true, the nonghost input fields have a smaller w-subscript than the
        # ghost input fields that they are brought into.  Therefore, e⁻ⁱᵏᴸ must be
        # multiplied to the nonghost input fields to create the ghost input fields.
        #
        # For isfwd = false, the nonghost input fields have a larger w-subscript than the
        # ghost input fields that they are brought into.  Therefore, e⁺ⁱᵏᴸ must be
        # multiplied to the nonghost input fields to create the ghost input fields.
        iw = isfwd ? Nw : 1
        Vₛ[Base.setindex(axes(Vₛ), iw:iw, nw)...] .*= e⁻ⁱᵏᴸ^ns
    else  # symmetry boundary
        # A. isfwd = true (forward difference)
        #
        # A-0. Goal
        # The forward difference is taken between primal grid planes.  Because the domain
        # boundaries are always primal grid planes, some input fields are on the domain
        # boundaries.  Those input fields are zero for the symmetry boundary condition and
        # do not contribute to the output derivatives.  Our goal is to make sure that the
        # forward derivative operator properly ignore those fields that should be zero.
        #
        # A-1. Determination of the indices of the the rows to modify
        # The negative-end boundary input fields contribute to the w-derivatives at the
        # first (nonghost) dual grid points in the w-direction.  The positive-end boundary
        # input fields contribute to the w-derivatives at the last dual grid points in the
        # w-direction.  Therefore, the array subscripts (i,j,k) we need to consider are the
        # ones with iw ≡ (i,j,k)[w] = 1 and Nw.
        #
        # A-2. Determination of the column in each row to modify
        # The negative-end boundary input fields have subscripts the same as the derivatives
        # they contribute to, so for iw = 1 we modify the diagonal entries stored in V₀.
        # The positive-end boundary input fields have w-subscript greater by 1 than the
        # derivatives they contribute to, so for iw = Nw we modify the off-diagonal entries
        # stored in Vₛ.
        #
        # A-3. Modification of the values in V
        # The negative-end boundary input fields must be zeroed, so for (i,j,k) with iw = 1
        # we set V₀[i,j,k] = 0.
        # The positive-end boundary input fields must be also zeroed, so for (i,j,k) with iw = Nw
        # we set Vₛ[i,j,k] = 0.  Note that the positive-end boundary input fields are ghost
        # fields that are brought from the negative-end boundary (nonghost) input fields.
        #
        #
        # B. isfwd = false (backward difference)
        #
        # B-0. Goal
        # The backward difference is taken between dual grid planes.  The domain boundaries
        # are always primal grid planes.  Therefore, no input fields are on the domain
        # boundaries, and thus no input fields needs to be zeroed.
        #
        # However, the fields at dual grid points are symmetric (rather than anti-symmetric)
        # around the symmetry boundaries (so that their interpolated values at the symmetry
        # boundary are nonzero).  Therefore, the difference between the two is zero, leading
        # to zero derivatives evaluated at the symmetry boundaries.
        #
        # On the positive-end boundary, the derivative of the dual field is not evaluated,
        # so we don't have to do anything to realize the positivte-end symmetry boundary.
        # On the other hand, on the negative-end boundary the darivative of the dual field
        # is evaluated.  We have to make sure that the resulting derivatives are zero.
        #
        # B-1. Determination of the indices of the the rows to modify
        # The output derivatives at the negative-end boundary are subscripted (i,j,k) with
        # iw ≡ (i,j,k)[w] = 1.
        #
        # B-2. Determination of the column in each row to modify
        # To make the output derivatives at the negative-end boundary zero, both input
        # fields contributing to the derivatives must be zeroed.  Therefore, both V₀ and Vₛ
        # must be modified.
        #
        # B-3. Modification of the values in V
        # We simply zero the values, i.e., V₀[i,j,k] = Vₛ[i,j,k] = 0, for (i,j,k) with iw = 1.
        #
        # The above procedure can be written compactly as follows without separating the
        # isfwd = true and false cases.

        # Zero the diagonal entries (V₀) multiplied with the fields on the boundary.
        #
        # Note that the code below is independent of isfwd.  When the grid is uniorm, the
        # operators for the same boundary condition but for the opposite isfwd must be the
        # minus of the transpose of each other, so the diagonals of two operators must have
        # zeros at the same locations.  This means the locations of zeros on the diagonal
        # must be independent of isfwd.
        iw = 1
        V₀[Base.setindex(axes(V₀), iw:iw, nw)...] .= 0

        # Zero the off-diagonal entries (Vₛ) multiplied with the fields on the boundary.
        #
        # Note that the code below is indepent of the boundary condition.  When the grid is
        # uniform, this means that for the same ns, the operators for the symmetry boundary
        # conditions will have the same off-diagonal entries, except for the opposite signs.
        iw = isfwd ? Nw : 1
        Vₛ[Base.setindex(axes(Vₛ), iw:iw, nw)...] .= 0
    end

    # Vectors with length M
    i = reshape(I, M)
    jₛ = reshape(Jₛ, M)
    v₀ = reshape(V₀, M)
    vₛ = reshape(Vₛ, M)

    # Vectors with length 2M
    Is = [i; i]  # row indices of [diagonal; off-diagonal] entries
    Js = [i; jₛ]  # column indices of [diagonal; off-diagonal] entries
    Vs = [v₀; vₛ]  # matrix entries of [diagonal, off-diagonal] entries

    return Is, Js, Vs
end


text:
# Pawel Zielinski
module matrixgen

using LinearAlgebra

export  blockmat

	function matcond(n::Int, c::Float64)
		# Function generates a random square matrix A of size n with
		# a given condition number c.
		# Inputs:
		#	n: size of matrix A, n>1
		#	c: condition of matrix A, c>= 1.0
		#
		# Usage: matcond(10, 100.0)
		#
		# Pawel Zielinski
        if n < 2
         error("size n should be > 1")
        end
        if c< 1.0
         error("condition number  c of a matrix  should be >= 1.0")
        end
        (U,S,V)=svd(rand(n,n))
        return U*diagm(0 =>[LinRange(1.0,c,n);])*V'
	end

  function blockmat(n::Int, l::Int, ck::Float64, outputfile::String)
		# Function generates a random block sparse matrix A of size n with
		# a given condition number ck of inner block Ak and it stores the output
		# matrix in a text file.
		# Inputs:
		#	n: size of block matrix A, n>3
		#   l: size of inner matrices Ak, n mod l =0 (n is  divisible by l), l>1
		#	ck: condition of inner matrix Ak, ck>= 1.0
		# outputfile: name of the output text file
		#
		# Usage: blockmat(100, 4 ,10.0, "A.txt")
		#
		#
		#  the output file format
	  #  n  l              <--- the size of block matrix A, the size of inner matrices Ak
		#  i1  j1   A[i1,j1] <--- a non-zero element of block matrix A
		#  i2  j2   A[i2,j2] <--- a non-zero element of block matrix A
		#  i3  j3   A[i3,j3] <--- a non-zero element of block matrix A
		#  ...
		#  ...
		#  EOF
		#

    if n < 4
     error("size n should be > 3")
    end
    if l < 2
     error("size l should be > 1")
    end
		if n%l!=0
			error("n is not divisible by l")
		end

		nb=div(n,l)
		Ak=Matrix{Float64}(undef, l, l)
		open(outputfile, "w") do f
			println(f, n," ",l)
			for k in 1:nb
				Ak=matcond(l, ck)
				for i in 1:l, j in 1:l
					println(f,(k-1)*l+i," ",(k-1)*l+j," ", Ak[i,j])
				end
				if k<nb
			  	 for i in 1:l
						println(f,(k-1)*l+i," ",k*l+i," ",0.3*rand())
				 	 end
				end
				if k>1
			   for i in 1:l
					 println(f,(k-1)*l+i," ",(k-1)*l-1," ",0.3*rand())
					 println(f,(k-1)*l+i," ",(k-1)*l," ",0.3*rand())
				 end
				end
			end
		end	 # do
	end # blockmat

end # matrixgen


text:
using TensorIntegration
using Test
using LinearAlgebra

@testset "Expectation of uniform by simpson rule, 1D/2D/3D" begin
    # add test for basic integrals
    N = 5
    # uniform expectation
    a = 1.
    b = 3.
    pdf(x, a, b) = 1.0 ./ (b-a)
    grid, weights = TensorIntegration.simpson(a, b, N)

    Exp = dot((grid.*pdf.(grid, a, b)),weights)
    @test Exp ≈ (a+b)/2.

    # now in 2d
    N2d = [5, 5]
    a2d = [1., 1.]
    b2d = [2., 2.]
    pdf(x, a, b, d) = 1.0 ./ (b[d]-a[d])

    grid, weights = TensorIntegration.tensor_simpson(a2d, b2d, N2d)
    Exp1 = dot((grid[:,1].*pdf(grid, a2d, b2d, 1)), weights)
    Exp2 = dot((grid[:,2].*pdf(grid, a2d, b2d, 2)), weights)

    @test Exp1[1] ≈ (a2d[1]+b2d[1])/2.
    @test Exp2[1] ≈ (a2d[2]+b2d[2])/2.

    # now in 3d
    N3d = [5, 5, 5]
    a3d = [1., 1., 0.1]
    b3d = [2., 2., 3.3]
    pdf(x, a::Array, b::Array) = 1.0 ./ prod(b3d.-a3d)

    grid, weights = TensorIntegration.tensor_simpson(a3d, b3d, N3d)
    Exp1 = dot((grid[:,1].*pdf(grid, a3d, b3d)), weights)
    Exp2 = dot((grid[:,2].*pdf(grid, a3d, b3d)), weights)
    Exp3 = dot((grid[:,3].*pdf(grid, a3d, b3d)), weights)

    @test isapprox(Exp1, (a3d[1]+b3d[1])/2.)
    @test isapprox(Exp2, (a3d[2]+b3d[2])/2.)
    @test isapprox(Exp3, (a3d[3]+b3d[3])/2.)

    # now with one dimension only 1 point
    N3d = [5, 5, 1]
    a3d = [1., 1., 1.]
    b3d = [2., 2., 2.]
    pdf(x, a::Array, b::Array) = 1.0 ./prod(b3d.-a3d)

    grid, weights = TensorIntegration.tensor_simpson(a3d, b3d, N3d)
    Exp1 = dot((grid[:,1].*pdf(grid, a3d, b3d)), weights)
    Exp2 = dot((grid[:,2].*pdf(grid, a3d, b3d)), weights)
    Exp3 = dot((grid[:,3].*pdf(grid, a3d, b3d)), weights)

    @test isapprox(Exp1, (a3d[1]+b3d[1])/2.)
    @test isapprox(Exp2, (a3d[2]+b3d[2])/2.)
    @test isapprox(Exp3, (a3d[3]+b3d[3])/2.)
end