text:
# module BlackBody

# include("Radiation.jl")

# export bbody, set_temp_kev, set_temp_k, set_temp_hz, set_lum, bb_spectrum
using Juriba.Constants
using Juriba.Radiation
const RD = Radiation

@Base.kwdef mutable struct BBody <: AbstractRadiation
    size::Int					        # size of arrays
    en_phot::Vector{Float64}			# array of photon energies
    num_phot::Vector{Float64} 			# array of number of photons in units of erg/s/Hz
    en_phot_obs::Vector{Float64}		# same as above but in observer frame
    num_phot_obs::Vector{Float64} 		# same as above but in observer frame

    z::Float64 = 0.				            # Dimensions of emitting region
    r::Float64 = 0.
    vol::Float64 = 0.					    # Volume of emitting region
    beta::Float64 = 0.				        # speed of the emitting region
    angle::Float64 = 0.		        # Viewing angle/Doppler factor of emitting region
    dopfac::Float64 = 0.
    dopnum::Float64 = 0.				        # Doppler boosting exponent, depends on geometry
    counterjet::Bool = false			# boolean switch if user wants to include counterjet emission
    geometry::String = ""   	        # string to track geometry of emitting region

    Tbb::Float64 = 0.
    Lbb::Float64 = 0.
    normbb::Float64 = 0.

end

function bbody(s=40)
    size = s

    en_phot = zeros(Float64, size)
    num_phot = zeros(Float64, size)
    en_phot_obs = zeros(Float64, size)
    num_phot_obs = zeros(Float64, size)

    return BBody(size=size, en_phot=en_phot, num_phot=num_phot, en_phot_obs=en_phot_obs, num_phot_obs=num_phot_obs)

end

function set_temp_kev!(bbody::BBody, T)

    Tbb = T*kboltz_kev2erg/kboltz

    emin = 0.02*kboltz*Tbb
    emax = 30.0*kboltz*Tbb

    einc = (log10(emax)-log10(emin))/(bbody.size-1)

    for i in eachindex(bbody.en_phot)
        bbody.en_phot[i] = 10^(log10(emin) + i*einc)
        bbody.en_phot_obs[i] = bbody.en_phot[i]
    end

    bbody.Tbb = Tbb
end

function set_temp_k!(bbody::BBody, T)

    Tbb = T

    emin = 0.02*kboltz*Tbb
    emax = 30.0*kboltz*Tbb

    einc = (log10(emax)-log10(emin))/(bbody.size-1)

    for i in eachindex(bbody.en_phot)
        bbody.en_phot[i] = 10^(log10(emin) + i*einc)
        bbody.en_phot_obs[i] = bbody.en_phot[i]
    end

    bbody.Tbb = Tbb
end

function set_temp_hz!(bbody::BBody, nu)
    Tbb = (herg*nu)/(2.82*kboltz)

    emin = 0.02*kboltz*Tbb
    emax = 30.0*kboltz*Tbb

    einc = (log10(emax)-log10(emin))/(bbody.size-1)

    for i in eachindex(bbody.en_phot)
        bbody.en_phot[i] = 10^(log10(emin) + i*einc)
        bbody.en_phot_obs[i] = bbody.en_phot[i]
    end

    bbody.Tbb = Tbb

end

function set_lum!(bbody::BBody, L)
    bbody.Lbb = L
    bbody.normbb = bbody.Lbb/(sbconst*bbody.Tbb^4)
end

function bb_spectrum!(bbody::BBody)
    for i ∈ 1:bbody.size
	    bbody.num_phot[i] = bbody.normbb*2π*herg*bbody.en_phot_obs[i]/herg^3.0/
                            (cee^2*(exp(bbody.en_phot_obs[i]/(bbody.Tbb*kboltz)) - 1.))
        bbody.num_phot_obs[i] = bbody.num_phot[i]
    end
end

function temp_kev(bbody::BBody)
    return bbody.Tbb*kboltz/kboltz_kev2erg
end

function temp_k(bbody::BBody)
    return bbody.Tbb
end

function temp_hz(bbody::BBody)
    return 2.82*kboltz*bbody.Tbb/herg
end

function lum(bbody::BBody)
    return bbody.Lbb
end

function norm(bbody::BBody)
    return bbody.normbb
end

function Urad(bbody::BBody, d)
    return bbody.Lbb/(4π*d^2*cee)
end

function test(bbody::BBody)
    println("Black body temperature in K: $(temp_k(bbody)) in keV: $(temp_kev(bbody))")
    println("Black body bolometric luminosity in erg s⁻¹: $(bbody.Lbb)")
end


text:
"""
  split_train_test(X::Array{Array{Float64,1},1}, y::Array; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)

Function to split X and Y in between a train and a test sample.
"""
function split_train_test(X::Array{Array{Float64,1},1}, y::Array{Array{Float64,1},1}; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)

    uni = Bernoulli(ratio)
    draws = rand(uni, size(X,1))

    XTest = X[draws .== 1]
    XTrain = X[draws .== 0]

    yTest = y[draws .== 1]
    yTrain = y[draws .== 0]

    if testRatio == true
        println("ratio XTest/XTrain = $(size(XTest,1) / size(XTrain,1))")
        println("ratio yTest/yTrain = $(size(yTest,1) / size(yTrain,1))")

        if abs(size(XTest,1)/size(XTrain,1) - ratio) > tolTestRatio
            Base.error("Splitting does not respect the split ratio for X.")
        end
        if abs(size(yTest,1)/size(yTrain,1) - ratio) > tolTestRatio
            Base.error("Splitting does not respect the split ratio for y.")
        end

    end

    return XTrain, XTest, yTrain, yTest

end

"""
  split_train_test(X::Array{Array{Float64,1},1}, y::Array{Array{Float64,1},1}, yConvergence::Array{Int64,1}; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)

Function to split X, y and yConvergence in between a train and a test sample.
"""
function split_train_test(X::Array{Array{Float64,1},1}, y::Array{Array{Float64,1},1}, yConvergence::Array{Int64,1}; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)

    uni = Bernoulli(ratio)
    draws = rand(uni, size(X,1))

    XTest = X[draws .== 1]
    XTrain = X[draws .== 0]

    yTest = y[draws .== 1]
    yTrain = y[draws .== 0]

    yConvergenceTest = yConvergence[draws .== 1]
    yConvergenceTrain = yConvergence[draws .== 0]

    if testRatio == true

        println("ratio XTest/XTrain = $(size(XTest,1) / size(XTrain,1))")
        println("ratio yTest/yTrain = $(size(yTest,1) / size(yTrain,1))")
        println("ratio yConvergenceTTest/yConvergenceTTrain = $(size(yConvergenceTTest,1) / size(yConvergenceTTrain,1))")

        if abs(size(XTest,1)/size(XTrain,1) - ratio) > tolTestRatio
            Base.error("Splitting does not respect the split ratio for X.")
        end
        if abs(size(yTest,1)/size(yTrain,1) - ratio) > tolTestRatio
            Base.error("Splitting does not respect the split ratio for y.")
        end
        if abs(size(yConvergenceTest,1)/size(yConvergenceTrain,1) - ratio) > tolTestRatio
            Base.error("Splitting does not respect the split ratio for yConvergence.")
        end

    end

    return XTrain, XTest, yTrain, yTest, yConvergenceTrain, yConvergenceTest

end

"""
  split_train_test(X::Array{Float64,2}, y::Array{Float64,2}; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)

Function to split X and Y in between a train and a test sample.
"""
function split_train_test(X::Array{Float64,2}, y::Array{Float64,2}; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)
  return split_train_test(convert_to_array_of_arrays(X), convert_to_array_of_arrays(y), ratio = ratio, testRatio = testRatio, tolTestRatio = tolTestRatio)
end

"""
  split_train_test(X::Array{Float64,2}, y::Array{Float64,2}, yConvergence::Array{Int64,1}; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)

Function to split X, y and yConvergence in between a train and a test sample.
"""
function split_train_test(X::Array{Float64,2}, y::Array{Float64,2}, yConvergence::Array{Int64,1}; ratio::Float64 = 0.25, testRatio::Bool=false, tolTestRatio::Float64 = 0.10)
  return return split_train_test(convert_to_array_of_arrays(X), convert_to_array_of_arrays(y), yConvergence, ratio = ratio, testRatio = testRatio, tolTestRatio = tolTestRatio)
end

"""
  convert_to_array_of_arrays(X::Array{Float64,2})

Function to convert a matrix of type X::Array{Float64,2} to an array of arrays
Array{Array{Float64,1},1}. The latter is the format expected by Scikit-Learn.jl
"""
function convert_to_array_of_arrays(X::Array{Float64,2})

  # initialization
  Y = Array{Array{Float64,1},1}(size(X,1))

  for i=1:size(X,1)
    Y[i] = X[i,:]
  end

  return Y

end

"""
  convert_to_array(X::Array{Array{Float64,1},1})

Function to convert an array of arrays Array{Array{Float64,1},1}
to a matrix of type X::Array{Float64,2} to.
"""
function convert_to_array(X::Array{Array{Float64,1},1})

  # initialization
  Y = zeros(size(X,1), size(X[1],1))

  for j=1:size(X[1],1)
  for i=1:size(X,1)
    Y[i,j] = X[i][j]
  end
  end

  return Y

end

"""
  convert_to_vector(X::Array{Array{Float64,1},1})

Function to convert an array of arrays Array{Array{Float64,1},1}
to a matrix of type X::Array{Float64,1} to.
"""
function convert_to_vector(X::Array{Array{Float64,1},1})

  #safety checks
  if size(X[1],1) > 1
    info("size(X,1) = $(size(X,1))")
    info("X[1] = $(X[1])")
    Base.error("size(X[1],1) = $(size(X[1],1)). Use convert_to_array.")
  end
  # initialization
  Y = zeros(size(X,1))

  for i=1:size(X,1)
    Y[i] = X[i][1]
  end

  return Y

end

"""
  convert_vector_to_array(X::Array{Float64,1})

Scikit-Learn expect the input of the predict function to be of type Array{Float64,2}.
This function converts an Array{Float64,1} to a one-observation Array{Float64,2}.
"""
function convert_vector_to_array(X::Array{Float64,1})

  Y = Array{Float64}(1, length(X))

  for i=1:length(X)
    Y[1,i] = X[i]
  end

  return Y

end

"""
  convert_array_to_vector(X::Array{Float64,2})

Function to convert a one-observation Array{Float64,2} to an Array{Float64,1}
"""
function convert_vector_to_array(X::Array{Float64,2})

  Y = Array{Float64,1}(size(X,2))

  for i=1:size(X,2)
    Y[i] = X[1,i]
  end

  return Y

end

"""
predict_scaled_robust(clfr::PyObject, clf::PyObject, XScaled::Array;
                              nonConvergenceFlag::Int64 = 0,
                              convergenceFlag::Int64 = 1,
                              penaltyValue::Float64 = 999.0)

Function to predict several observations. Robust in the sense that the classifier
is first use to determine whether or not the model converges.
* If the model converges: value is given by the regression model, trained on data
  for which convergence was reached
* If the model does not converge: value = penalty value
"""
function predict_scaled_robust(clfr::PyObject, clf::PyObject, XScaled::Array;
                              nonConvergenceFlag::Int64 = 0,
                              convergenceFlag::Int64 = 1,
                              penaltyValue::Float64 = 999.0)

  # Classifier:
  #------------
  convergenceChek = round(Int64, predict(clf, XScaled))

  # Regression model:
  #------------------
  predictedValues = predict(clfr, XScaled)

  # When the model is not expected to converge, set value = penalty value
  #----------------------------------------------------------------------
  for i=1:size(predictedValues,1)
    if convergenceChek[i] == nonConvergenceFlag
      predictedValues[i, :] = penaltyValue
    end
  end

  return predictedValues

end

"""
  predict_one_obs_scaled(clfr::PyObject, XScaled::Array{Float64,1})

Function to predict one observation. Assume that X has been scaled accordingly.
"""
function predict_one_obs_scaled(clfr::PyObject, XScaled::Array{Float64,1})

  convert_vector_to_array(predict(clfr, convert_vector_to_array(XScaled)))

end

"""
  predict_one_obs_unscaled(clfr::PyObject, XScaled::Array{Float64,1})

Function to predict one observation. Scale the observation first.
"""
function predict_one_obs_unscaled(clfr::PyObject, X::Array{Float64,1}, scaler::PyObject)

  XScaled = transform(scaler, convert_vector_to_array(X))
  convert_vector_to_array(predict(clfr, convert_vector_to_array(XScaled)))

end

"""
  predict_one_obs_unscaled(clfr::PyObject, X::Array{Array{Float64,1}}, scaler::PyObject)

Function to predict one observation. Scale the observation first.
"""
function predict_one_obs_unscaled(clfr::PyObject, X::Array{Array{Float64,1}}, scaler::PyObject)

  XScaled = transform(scaler, X)
  convert_vector_to_array(predict(clfr, XScaled))

end

"""
  set_model_function!(sModelsProblem::SModelsProblem, f::Function)

Function to set the field sModelsProblem.modelFunction
"""
function set_model_function!(sModelsProblem::SModelsProblem, f::Function)

  sModelsProblem.modelFunction = f

end

"""
  set_bounds!(sModelsProblem::SModelsProblem, lowerBound::Array{Float64,1}, upperBound::Array{Float64,1})
"""
function set_bounds!(sModelsProblem::SModelsProblem, lowerBound::Array{Float64,1}, upperBound::Array{Float64,1})

  sModelsProblem.lowerBound = lowerBound
  sModelsProblem.upperBound = upperBound

end

"""
  calculate_mean_per_error(yTrue::Array{Float64,1}, yPredicted::Array{Float64,1})

Function to calculate the mean percentage error
"""
function calculate_mean_per_error(yTrue::Array{Float64,1}, yPredicted::Array{Float64,1})

  perError = zeros(length(yTrue))
  meanPredicted = mean(yPredicted)

  for i=1:length(yTrue)
    if yPredicted[i] != 0
      perError[i] = (yPredicted[i] - yTrue[i])/yTrue[i]
    else
      perError[i] = (yPredicted[i] - yTrue[i])/meanPredicted
    end
  end

  return mean(perError)

end

"""
  calculate_mean_per_error(yTrue::Array{Float64,2}, yPredicted::Array{Float64,2})

Function to calculate the mean percentage error
"""
function calculate_mean_per_error(yTrue::Array{Float64,2}, yPredicted::Array{Float64,2})

  perError = zeros(size(yTrue, 1),size(yTrue, 2))
  meanPredicted = mean(yPredicted)

  for j=1:size(yTrue, 2)
    for i=1:size(yTrue, 1)
      if yPredicted[i,j] != 0
        perError[i,j] = (yPredicted[i,j] - yTrue[i,j])/yTrue[i,j]
      else
        perError[i,j] = (yPredicted[i,j] - yTrue[i,j])/meanPredicted[j]
      end

    end
  end

  return mean(mean(perError))

end

"""
  calculate_mean_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,2})

Function to calculate the mean percentage error
"""
function calculate_mean_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,2})

  return calculate_mean_per_error(convert_to_array(yTrue), yPredicted)

end

"""
  calculate_mean_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,1})

Function to calculate the mean percentage error
"""
function calculate_mean_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,1})

  return calculate_mean_per_error(convert_to_vector(yTrue), yPredicted)

end

"""
  calculate_maximum_abs_per_error(yTrue::Array{Float64,1}, yPredicted::Array{Float64,1})

Function to calculate the maximum of the absolute value fo the percentage difference
"""
function calculate_maximum_abs_per_error(yTrue::Array{Float64,1}, yPredicted::Array{Float64,1})

  perError = zeros(length(yTrue))
  meanPredicted = mean(yPredicted)

    for i=1:size(yTrue, 1)
      if yPredicted[i] != 0
        perError[i] = abs((yPredicted[i] - yTrue[i])/yTrue[i])
      else
        perError[i] = abs((yPredicted[i] - yTrue[i])/meanPredicted)
      end
    end

  return maximum(perError)

end

"""
  calculate_maximum_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,1})

Function to calculate the maximum of the absolute value fo the percentage difference
"""
function calculate_maximum_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,1})

  return calculate_maximum_abs_per_error(convert_to_vector(yTrue), yPredicted::Array{Float64,1})

end

"""
  calculate_maximum_abs_per_error(yTrue::Array{Float64,2}, yPredicted::Array{Float64,2})

Function to calculate the maximum of the absolute value fo the percentage difference
"""
function calculate_maximum_abs_per_error(yTrue::Array{Float64,2}, yPredicted::Array{Float64,2})

  perError = zeros(size(yTrue, 1),size(yTrue, 2))
  meanPredicted = mean(yPredicted)

  for j=1:size(yTrue, 2)
    for i=1:size(yTrue, 1)
      if yPredicted[i,j] != 0
        perError[i,j] = abs((yPredicted[i,j] - yTrue[i,j])/yTrue[i,j])
      else
        perError[i,j] = abs((yPredicted[i,j] - yTrue[i,j])/meanPredicted[j])
      end

    end
  end

  return maximum(perError)

end

"""
  calculate_maximum_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,2})

Function to calculate the maximum of the absolute value fo the percentage difference
"""
function calculate_maximum_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,2})

  return calculate_maximum_abs_per_error(convert_to_array(yTrue), yPredicted)

end

"""
  calculate_median_abs_per_error(yTrue::Array{Float64,2}, yPredicted::Array{Float64,2})

Function to calculate the median of the absolute value fo the percentage difference
"""
function calculate_median_abs_per_error(yTrue::Array{Float64,2}, yPredicted::Array{Float64,2})

  perError = zeros(size(yTrue, 1),size(yTrue, 2))
  meanPredicted = mean(yPredicted)

  for j=1:size(yTrue, 2)
    for i=1:size(yTrue, 1)
      if yPredicted[i,j] != 0
        perError[i,j] = abs((yPredicted[i,j] - yTrue[i,j])/yTrue[i,j])
      else
        perError[i,j] = abs((yPredicted[i,j] - yTrue[i,j])/meanPredicted[j])
      end

    end
  end

  return median(perError)

end

"""
  calculate_median_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,2})

Function to calculate the median of the absolute value fo the percentage difference
"""
function calculate_median_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,2})

  return calculate_median_abs_per_error(convert_to_array(yTrue), yPredicted)

end

"""
  calculate_median_abs_per_error(yTrue::Array{Float64,1}, yPredicted::Array{Float64,1})

Function to calculate the median of the absolute value fo the percentage difference
"""
function calculate_median_abs_per_error(yTrue::Array{Float64,1}, yPredicted::Array{Float64,1})

  perError = zeros(length(yTrue))
  meanPredicted = mean(yPredicted)

    for i=1:length(yTrue)
      if yPredicted[i] != 0
        perError[i] = abs((yPredicted[i] - yTrue[i])/yTrue[i])
      else
        perError[i] = abs((yPredicted[i] - yTrue[i])/meanPredicted)
      end

    end

  return median(perError)

end

"""
  calculate_median_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,1})

Function to calculate the median of the absolute value fo the percentage difference
"""
function calculate_median_abs_per_error(yTrue::Array{Array{Float64,1},1}, yPredicted::Array{Float64,1})

  return calculate_median_abs_per_error(convert_to_vector(yTrue), yPredicted)

end

"""
  date_now()

Return date and time
"""
function date_now()

  "$(Dates.today())--$(Dates.hour(Dates.now()))h-$(Dates.minute(Dates.now()))m-$(Dates.second(Dates.now()))s"

end

ext: jl


text:
using ViewerGL, LinearAlgebraicRepresentation
GL = ViewerGL
Lar = LinearAlgebraicRepresentation

# setup structure
V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);
cuboid_big = (V,EV,FV,CV);
V,(VV,EV,FV,CV) = Lar.cuboid([0.75,0.75,0.75],true,[0.25,0.25,0.25]);
cuboid_small = (V,EV,FV,CV);

str = Lar.Struct([ cuboid_big, cuboid_small ]);
V,EV,FV,CV = Lar.struct2lar(str);

GL.VIEW([ GL.GLGrid(V,EV, GL.COLORS[1]), GL.GLFrame ]);
GL.VIEW([ GL.GLPol(V,CV, GL.COLORS[1]), GL.GLFrame ]);

# V, (copEV, copFE, copCF), boolmatrix = Lar.bool3d(str)

# space arrangement
cop_EV = Lar.coboundary_0(EV);
cop_FE = Lar.coboundary_1(V, FV, EV);
W = convert(Lar.Points, V');
W, copEV, copFE, copCF = Lar.spatial_arrangement(W, cop_EV, cop_FE);

# compute containment graph of components
bicon_comps = Lar.Arrangement.biconnected_components(copEV)

# compute euler characteristic
χ = Lar.euler_characteristic(W, copEV, copFE)

# triangulation
triangulated_faces = Lar.triangulate(W, [copEV, copFE]);
FVs = convert(Array{Lar.Cells}, triangulated_faces);
GL.VIEW( GL.GLExplode(V,FVs, 2.,2.,2.,99,1) );

# v,cvs,fvs,evs = Lar.pols2tria(V, copEV, copFE, copCF) # V by cols
# GL.VIEW( GL.GLExplode(V,cvs[2:end],1.5,1.5,1.5,99) );

EVs = Lar.FV2EVs(copEV, copFE); # polygonal face fragments
GL.VIEW( GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1) );


text:
using DFControl, Test

testjobpath = joinpath(testdir, "testassets", "test_job")
job = DFJob(testjobpath);

#output stuff
out = outputdata(job;print=false,onlynew=false);
@test haskey(out, "nscf")
@test haskey(out["nscf"], :fermi)

@test isconverged(job["scf"])
nscf = DFControl.input(job, "nscf")
nscf2 = DFInput(nscf, "nscf2", data=[:testdata => (:testoption, "test"), :k_points => (:blabla, [1,1,1,1,1,1])])

@test data(nscf2, :testdata).option == :testoption
@test data(nscf2, :testdata).data   == "test"
@test data(nscf2, :testdata).name   == :testdata
@test data(nscf2, :k_points).option == :blabla
@test data(nscf2, :k_points).data   == [1,1,1,1,1,1]
@test data(job["scf"], :k_points).option == :automatic
@test_throws ErrorException job["bladbsflk"]
@test_throws ErrorException nscf2[:bladkfj]

setkpoints!(nscf2, (3,3,3), print=false)
@test data(nscf2, :k_points).data  == DFControl.kgrid(3, 3, 3, :nscf)

setkpoints!(nscf2, [(3.,3.,3.,1.), (3.,3.,3.,1.)], print=false)
@test data(nscf2, :k_points).option  == :crystal_b
@test data(nscf2, :k_points).data  == [(3.,3.,3.,1.), (3.,3.,3.,1.)]

setkpoints!(nscf2, (3,3,3,0,0,1), print=false)
@test data(nscf2, :k_points).option  == :automatic
@test data(nscf2, :k_points).data  == [3,3,3,0,0,1]

fermi = outputdata(job, "nscf")[:fermi]
@test fermi == 17.4572

push!(job, gencalc_bands(job["bands"], [(0.5,0.0,0.5,10.0),(0.0,0.0,0.0,10.0),(0.5,0.5,0.5,1.0)], name="bands2"))
@test job["bands2"][:calculation] == "bands"
@test data(job, "bands2", :k_points).data == [(0.5,0.0,0.5,10.0),(0.0,0.0,0.0,10.0),(0.5,0.5,0.5,1.0)]
push!(job, gencalc_nscf(job["nscf"], (10,10,10), name="nscf2"))
@test job["nscf2"][:calculation] == "nscf"
push!(job, gencalc_scf(job["scf"], (5,5,5,1,1,1), name="1scf2"))
@test job["1scf2"][:calculation] == "scf"

wanflags = [:write_hr => true, :wannier_plot => true]

setprojections!(job, :Pt => [:s, :p, :d])
wancalc = gencalc_wan(job["nscf"], structure(job), fermi-7.0; Epad=5.0, wanflags=wanflags)
push!.((job,), wancalc)
@test job["wan"][:write_hr] == job["wan"][:wannier_plot] == true

wanout = outputdata(job, "wan")
@test length(wanout[:final_state]) == 20
@test length(wanout[:wannierise]) == 203
job.inputs = job.inputs[1:4]

setflags!(job, :nspin => 2, print=false)
@test job["nscf"][:nspin] == 2
job["nscf"][:nspin] = 3
@test job["nscf"][:nspin] == 3

job[:nspin] = 2
@test job["nscf"][:nspin] == 2

job["nscf"][:Hubbard_U] = [1.0]
@test job["nscf"][:Hubbard_U] == [1.0]

job["nscf"][:starting_magnetization] = [[1.0, 2.0, 1.0]]
@test job["nscf"][:starting_magnetization] == [[1.0, 2.0, 1.0]]

job["nscf"][:Hubbard_J] = [0 1 2]
@test job["nscf"][:Hubbard_J] == [0 1 2]
push!.((job,), gencalc_wan(nscf, structure(job), fermi-7.0, Epad=5.0, wanflags=wanflags))

#TODO: add test with the new wancalc from projwfc

setflow!(job, ""=>false)
@test job.inputs[1].run == false
setflow!(job, "nscf" => true, "bands" => true)
@test job.inputs[3].run

@show job.inputs
save(job)
job2 = DFJob(local_dir)
@show job2.inputs
begin
    for (calc, calc2) in zip(job.inputs, job2.inputs)

        for (f, v) in calc.flags
            if f in (:Hubbard_J, :pseudo_dir)
                continue
            else
                @test v == calc2.flags[f]
            end
        end
        for (b1, b2) in zip(calc.data, calc2.data)
            for n in fieldnames(typeof(b1))
                @test getfield(b1, n) == getfield(b2, n)
            end
        end
    end
end

job3 = DFJob(job2, :lspinorb => true)
@test all(atoms(job3).==atoms(job2))
@test length(job3[:lspinorb]) == length(searchinputs(job3, QE))
rmflags!(job3, :lspinorb, print=false)

begin
    for (calc, calc2) in zip(job.inputs, job3.inputs)
        for (f, v) in calc.flags
            if f in (:Hubbard_J, :pseudo_dir)
                continue
            else
                @test v == calc2.flags[f]
            end
        end
        for (b1, b2) in zip(calc.data, calc2.data)
            for n in fieldnames(typeof(b1))
                @test getfield(b1, n) == getfield(b2,n)
            end
        end
    end
end

setcutoffs!(job)
@test job["scf"][:ecutwfc] == 32.0

setpseudos!(job, :Pt => Pseudo("Pt.UPF", "testassets/pseudos"))
@test job.structure.atoms[1].pseudo == Pseudo("Pt.UPF", "testassets/pseudos")
setpseudos!(job, :Pt, :test)
@test job.structure.atoms[1].pseudo == Pseudo("Pt.UPF", "testassets/pseudos")

setpseudos!(job, :test)
@test job.structure.atoms[1].pseudo == Pseudo("Pt.UPF", "testassets/pseudos")

testorbs = [:s, :p]
setprojections!(job, :Pt => testorbs)
@test convert.(Symbol, [p.orb for p in projections(job, :Pt)]) == testorbs
setwanenergies!(job, nscf, fermi-7.0, Epad=3.0)

@test job["wanup"][:dis_froz_max] == 13.2921
@test job["wanup"][:dis_win_max] == 16.292099999999998

setexecflags!(job, "pw.x", :nk => 230, :ndiag => 3)
@test DFControl.getfirst(x->x.symbol==:nk, execs(job, "nscf")[2].flags).value == 230
@test DFControl.getfirst(x->x.symbol==:ndiag, execs(job, "nscf")[2].flags).value == 3
rmexecflags!(job, "pw.x", :nk)
@test isempty(filter(x->x.symbol == :nk, execs(job, "nscf")[2].flags))

setexecdir!(job, "pw.x", joinpath(homedir(), "bin"))
@test execs(job, "nscf")[2].dir == joinpath(homedir(), "bin")

setname!(job, "nscf", "test")
@test DFControl.inpath(job, "test") == joinpath(job.local_dir, "test.in")
setname!(job, "test", "nscf")

setserverdir!(job, "localhost")
@test job.server_dir == "localhost"

setheaderword!(job, "defpart" => "frontend", print=false)
@test any(occursin.("frontend",job.header))

setdataoption!(job, "nscf",:k_points, :blabla, print=false)
@test data(job, "nscf", :k_points).option == :blabla
setdataoption!(job, :k_points, :test, print=false)
@test data(job, "nscf", :k_points).option == :test

report = progressreport(job; onlynew=false, print=false)
@test report[:fermi] == 17.4572
@test length(report[:accuracy]) == 9
oldat = atoms(job)[1]
newatompos = outputdata(job, "vc_relax", onlynew=false)[:final_structure]
job.structure = newatompos
push!(job.structure.atoms, oldat)

cp(joinpath(testdir, "testassets", "pseudos"), joinpath(testdir, "testassets", "pseudos_copy"), force=true)
setpseudos!(job, :Si => Pseudo("Si.UPF", joinpath(testdir, "testassets", "pseudos_copy")))
save(job)
@test ispath(joinpath(job.local_dir, "Si.UPF"))
@test atom(job, :Si).pseudo == Pseudo("Si.UPF", job.local_dir)
rm(joinpath(testdir, "testassets", "pseudos_copy"), recursive=true)

job["nscf"][:occupations] = "smearing"
job["nscf"][:degauss] = 2.0
job["nscf"][:smearing] = "mp"
projwfc = gencalc_projwfc(job["nscf"], -20, 10, 0.05)
@test projwfc[:Emin] == -20
@test projwfc[:degauss] == job["nscf"][:degauss]
@test projwfc[:ngauss] == 1

set_Hubbard_U!(job, :Si => 1.0)
@test atoms(job, :Si)[1].dftu.U == 1.0

prev_a = cell(job)[1, :]
prev_b = cell(job)[2, :]
prev_c = cell(job)[3, :]
prev_pos = position_cart.(atoms(job))
scale_cell!(job, 2)
@test prev_a .* 2 == cell(job)[1, :]
@test prev_b .* 2 == cell(job)[2, :]
@test prev_c .* 2 == cell(job)[3, :]
for (p, at) in zip(prev_pos, atoms(job))
	@test round.(DFControl.ustrip(p * 2), digits=5) == round.(DFControl.ustrip(position_cart(at)), digits=5)
end

set_magnetization!(job, :Pt => [1.0, 0.0, 0.0])
@test magnetization(atom(job, :Pt)) == DFControl.Vec3(1.0, 0.0, 0.0)

at = atoms(job)[1]
c = cell(job)
orig_pos_cart  = position_cart(at)
orig_pos_cryst = position_cryst(at)

set_position!(at, orig_pos_cart .+ c * [0.1, 0.1, 0.1], c)
@test isapprox(position_cryst(at), orig_pos_cryst .+ [0.1, 0.1, 0.1])

set_position!(at, orig_pos_cryst .+ [0.1, 0.1, 0.1], c)
@test position_cart(at) == c * position_cryst(at)

at2 = atoms(job)[2]
p1, p2 = position_cart(at), position_cart(at2)
mid = (p1+p2)/2
bondlength = distance(at, at2)
scale_bondlength!(at, at2, 0.5, c)

@test isapprox((position_cart(at) + position_cart(at2))/2, mid)
@test isapprox(distance(at, at2), bondlength/2)

@test isapprox(bandgap(job), 3.6522999999999994)

t = job["nscf"]
curlen = length(job.inputs)
id = findfirst(x->x.name == "nscf", job.inputs)
n = pop!(job, "nscf")
@test n == t
@test length(job.inputs) == curlen - 1

nscf = job["scf"]
rmflags!(job, :nspin, :lda_plus_u, :noncolin)
set_magnetization!(job, :Pt => [0.2, 1.0, 0.2])
DFControl.sanitizeflags!(job)
set_Hubbard_U!(job, :Pt => 2.3)
@test job["scf"][:lda_plus_u]
@test job["scf"][:noncolin]

rmflags!(job, :nspin, :lda_plus_u, :noncolin)
set_magnetization!(job, :Pt => [0.0, 0.0, 0.5])
DFControl.sanitizeflags!(job)
@test job["scf"][:nspin] == 2

using LinearAlgebra
new_str = create_supercell(structure(job), 1, 0, 0)
prevcell = cell(job)
@test norm(cell(new_str)[:,1]) == norm(prevcell[:, 1]) * 2
prevlen = length(atoms(job))
@test length(atoms(new_str)) == 2*prevlen
prevlen_Pt = length(atoms(job, :Pt))

set_magnetization!(job, :Pt => [0, 0, 1])
orig_projs = projections(atoms(job, :Pt)[1])

job.structure = create_supercell(structure(job), 1, 0, 0, make_afm=true)

DFControl.sanitize_magnetization!(job)
@test length(atoms(job, :Pt1)) == prevlen_Pt
@test length(atoms(job, :Pt2)) == prevlen_Pt

@test magnetization(atoms(job, :Pt1)[1]) == [0, 0, 1]
@test magnetization(atoms(job, :Pt2)[1]) == [0, 0, -1]

DFControl.sanitize_projections!(job)
@test projections(atoms(job, :Pt1)[1]) != projections(atoms(job, :Pt2)[1])

rm.(DFControl.inpath.(job.inputs))

rm(joinpath(splitdir(DFControl.inpath(job.inputs[1]))[1], "pw2wan_wanup.in"))
rm(joinpath(splitdir(DFControl.inpath(job.inputs[1]))[1], "pw2wan_wandn.in"))
rm(joinpath(job.local_dir, "job.tt"))
rm.(joinpath.((job.local_dir,), filter(x -> occursin("UPF", x), readdir(job.local_dir))))


text:
# Plot the various c_n components:

z = linspace(0,1,1000)

using PyPlot

clf()
plot(z,ones(z),label="n=0")

plot(z,z,label="n=1")

plot(z,4.*z.^2-2.,label="n=2")
plot(z,5.*z.^3-3.*z,label="n=3")
plot(z,6.*z.^4-4.*z.^2,label="n=4")
plot(z,7.*z.^5-5.*z.^3,label="n=5")
legend(loc="upper left")
savefig("plot_components.pdf", bbox_inches="tight")

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
# Potential energy function
function upair(r,data)
  r6 = r^6
  r12 = r6^2
  return data.eps4*(data.sig12/r12 - data.sig6/r6)
end


text:
using Combinatorics

const stirlings2_cache = Dict{Pair{Int, Int}, BigInt}()

function stirlings2(n::Int, k::Int)::BigInt
    if haskey(stirlings2_cache, Pair(n, k))
        return stirlings2_cache[Pair(n, k)]
    elseif n < 0
        throw(DomainError(n, "n must be nonnegative"))
    elseif n == k == 0
        return one(n)
    elseif n == 0 || k == 0
        return zero(n)
    elseif k == n - 1
        return binomial(n, 2)
    elseif k == 2
        return 2^(n - 1) - 1
    end

    result = k * stirlings2(n - 1, k) + stirlings2(n - 1, k - 1)
    stirlings2_cache[Pair(n, k)] = result
    return result
end

const factorial_cache = Dict{Int, BigInt}()

function cached_factorial(n::Int)::BigInt
    if haskey(factorial_cache, n)
        return factorial_cache[n]
    end

    result = factorial(BigInt(n))
    factorial_cache[n] = result
    return result
end

const fubini_cache = Dict{Int, BigInt}()

function fubini(n::Int)::BigInt
    if haskey(fubini_cache, n)
        return fubini_cache[n]
    end

    result = sum(cached_factorial(k) * stirlings2(n, k) for k in 0:n)
    fubini_cache[n] = result
    return result
end

function fubinis(n::Int)
    Channel{BigInt}(n) do ch
        for i in 1:n
            put!(ch, fubini(i))
        end
    end
end

function print_fubinis(n::Int)
    # @time begin
        print("[")
        for i in 1:n
            if i > 1
                print(", ")
            end
            print(fubini(i))
        end
        print("]")
        println()
    # end
end

if length(ARGS) < 1
    throw(ErrorException("Missing argument."))
end
n = parse(Int, ARGS[1])
print_fubinis(n)


text:
# Latafat, Themelis, Patrinos, "Block-coordinate and incremental aggregated
# proximal gradient methods for nonsmooth nonconvex problems."
# arXiv:1906.10053 (2019).
#
# Latafat. "Distributed proximal algorithms for large-scale structured optimization"
# PhD thesis, KU Leuven, 7 2020.
#
# Mairal, "Incremental majorization-minimization optimization with application to
# large-scale machine learning."
# SIAM Journal on Optimization 25, 2 (2015), 829–855.
#
# Defazio, Domke, "Finito: A faster, permutable incremental gradient method
# for big data problems."
# In International Conference on Machine Learning (2014), pp. 1125-1133.
#

using LinearAlgebra
using ProximalOperators
using ProximalAlgorithms.IterationTools
using Printf
using Base.Iterators
using Random
using StatsBase: sample

export solution

include("Finito_basic.jl")
include("Finito_LFinito.jl")
include("Finito_adaptive.jl")

struct Finito{R<:Real}
    γ::Maybe{Union{Array{R},R}}
    sweeping::Int8
    LFinito::Bool
    adaptive::Bool
    minibatch::Tuple{Bool,Int}
    maxit::Int
    verbose::Bool
    freq::Int
    α::R
    tol::R
    tol_b::R
    function Finito{R}(;
        γ::Maybe{Union{Array{R},R}} = nothing,
        sweeping = 1,
        LFinito::Bool = false,
        adaptive::Bool = false,
        minibatch::Tuple{Bool,Int} = (false, 1),
        maxit::Int = 10000,
        verbose::Bool = false,
        freq::Int = 10000,
        α::R = R(0.999),
        tol::R = R(1e-8),
        tol_b::R = R(1e-9),
    ) where {R}
        @assert γ === nothing || minimum(γ) > 0
        @assert maxit > 0
        @assert tol > 0
        @assert tol_b > 0
        @assert freq > 0
        new(γ, sweeping, LFinito, adaptive, minibatch, maxit, verbose, freq, α, tol, tol_b)
    end
end

function (solver::Finito{R})(
    x0::AbstractArray{C};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    N = N,
) where {R,C<:RealOrComplex{R}}

    stop(state) = false

    disp(it, state) = @printf "%5d | %.3e  \n" it state.hat_γ

    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    # dispatching the iterator
    if solver.LFinito
        iter = FINITO_LFinito_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
        )
    elseif solver.adaptive
        iter = FINITO_adaptive_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.tol,
            solver.tol_b,
            solver.sweeping,
            solver.α,
        )
    else
        iter = FINITO_basic_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
        )
    end

    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)

    num_iters, state_final = nothing, nothing
    for (it_, state_) in iter  # unrolling the iterator
        # see https://docs.julialang.org/en/v1/manual/interfaces/index.html
        if solver.verbose && mod(it_, solver.freq) == 0
            disp(it_, state_)
        end
        num_iters, state_final = it_, state_
    end
    if solver.verbose && mod(num_iters, solver.freq) !== 0
        disp(num_iters, state_final)
    end # for the final iteration
    return solution(state_final), num_iters
end

"""
    Finito([γ, sweeping, LFinito, adaptive, minibatch, maxit, verbose, freq, tol, tol_b])

Instantiate the Finito algorithm for solving fully nonconvex optimization problems of the form

    minimize 1/N sum_{i=1}^N f_i(x) + g(x)

where `f_i` are smooth and `g` is possibly nonsmooth, all of which may be nonconvex.

If `solver = Finito(args...)`, then the above problem is solved with

	solver(x0, [F, g, N, L])

where F is an array containing f_i's, x0 is the initial point, and L is an array of
smoothness moduli of f_i's; it is optional in the adaptive mode or if γ is provided.

Optional keyword arguments are:
* `γ`: an array of N stepsizes for each coordinate
* `sweeping::Int` 1 for uniform randomized (default), 2 for cyclic, 3 for shuffled
* `LFinito::Bool` low memory variant of the Finito/MISO algorithm
* `adaptive::Bool` to activate adaptive smoothness parameter computation
* `minibatch::(Bool,Int)` to use batchs of a given size
* `maxit::Integer` (default: `10000`), maximum number of iterations to perform.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `10000`), frequency of verbosity.
* `α::R` parameter where γ_i = αN/L_i
* `tol::Real` (default: `1e-8`), absolute tolerance for the adaptive case
* `tol_b::R` tolerance for the backtrack (default: `1e-9`)
"""

Finito(::Type{R}; kwargs...) where {R} = Finito{R}(; kwargs...)
Finito(; kwargs...) = Finito(Float64; kwargs...)

"""
If `solver = Finito(args...)`, then

    itr = iterator(solver, x0, [F, g, N, L])

is an iterable object. Note that [maxit, verbose, freq] fields of the solver are ignored here.

The solution at any given state can be obtained using solution(state), e.g.,
for state in Iterators.take(itr, maxit)
    # do something using solution(state)
end

See https://docs.julialang.org/en/v1/manual/interfaces/index.html
and https://docs.julialang.org/en/v1/base/iterators/ for a list of iteration utilities
"""

function iterator(
    solver::Finito{R},
    x0::AbstractArray{C};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    N = N,
) where {R,C<:RealOrComplex{R}}
    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    # dispatching the iterator
    if solver.LFinito
        iter = FINITO_LFinito_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
        )
    elseif solver.adaptive
        iter = FINITO_adaptive_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.tol,
            solver.tol_b,
            solver.sweeping,
            solver.α,
        )
    else
        iter = FINITO_basic_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
        )
    end
    return iter
end

max_issues_repo_name: kul-forbes/CIAOAlgorithms

max_forks_repo_name: kul-forbes/CIAOAlgorithms


text:
# read dial-a-ride problem files

#function to read vehicle data
function readVehicleFile(filename::String)

	tables = readTablesFromFile(filename)

	table = tables["vehicles"]
	n = size(table.data,1) # number of vehicles
	assert(n >= 1)

	# create vehicles from data in table
	vehicles = Vector{Vehicle}(n)
	columns = table.columns # shorthand
	for i = 1:n
		vehicles[i] = Vehicle()
		vehicles[i].index = columns["index"][i]
		vehicles[i].capacity = columns["capacity"][i]
		vehicles[i].cost_coeff = columns["cost_coeff"][i]
		vehicles[i].max_time = columns["max_time"][i]
		vehicles[i].speed = columns["speed"][i]
		assert(vehicles[i].index == i)
	end

	return vehicles
end

#function to read node data
function readNodeFile(filename::String)

	tables = readTablesFromFile(filename)

	table = tables["nodes"]
	n = size(table.data,1) #number of nodes
	assert(n >= 1)

	#create nodes from data in table
	nodes = Vector{Node}(n)
	columns = table.columns #shorthand
	for i = 1:n
		nodes[i] = Node()
		nodes[i].index = columns["index"][i]
		nodes[i].name = columns["name"][i]
		nodes[i].x_loc = columns["x_location"][i]
		nodes[i].y_loc = columns["y_location"][i]
	end

	return nodes

end

#function to read request data
function readRequestFile(filename::String)

	tables = readTablesFromFile(filename)

	table = tables["requests"]
	n = size(table.data,1) #number of requests
	assert(n >= 1)

	#create requests from data in table
	requests = Vector{Request}(n)
	columns = table.columns #shorthand
	for i = 1:n
		requests[i] = Request()
		requests[i].index = columns["index"][i]
		requests[i].num_pass = columns["num_pass"][i]
		requests[i].origin = Location(columns["origin_x_loc"][i],columns["origin_y_loc"][i])
		requests[i].destination = Location(columns["destination_x_loc"][i],columns["destination_y_loc"][i])
		#requests[i].from_x_loc = columns["from_x_loc"][i]
		#requests[i].from_y_loc = columns["from_y_loc"][i]
		#requests[i].to_x_loc = columns["to_x_loc"][i]
		#requests[i].to_y_loc = columns["to_y_loc"][i]
		requests[i].early = columns["early"][i]
		requests[i].late = columns["late"][i]
	end

	return requests

end

max_stars_repo_name: jprendo/TaaSSim


text:
export load_img, save_img, split_rgb

using Images

"""
    load_img(filename::String)

Load an image given file name.
"""
function load_img(filename::String)
    img = FileIO.load(filename)
    convert(Array{UInt8}, rawview(channelview(img)))
end

"""
    load_img(filename::String)

Load a pre-saved Atari image given ROM name and frame number.
"""
function load_img(rom_name::String, frame_number::Int64)
    filename = string(@__DIR__, "/../images/", rom_name, "_frame_30.png")
    load_img(filename)
end

"""
    load_rgb(rom_name::String, frame_number::Int64)

Load a pre-saved Atari image given ROM name and frame number.
Return an array containing the separated chanels.
"""
function load_rgb(rom_name::String, frame_number::Int64)
    filename = string(@__DIR__, "/../images/", rom_name, "_frame_$frame_number.png")
    split_rgb(load_img(filename))
end

"""
    save_img(img::Array{Float64,2}, filename::String)

Save an image to the input path.
"""
function save_img(img::Array{Float64}, filename::String)
    n_dims = ndims(img)
    if n_dims == 2
        FileIO.save(filename, colorview(Gray, img))
    elseif n_dims == 3
        FileIO.save(filename, colorview(RGB, img))
    else
        println("Warning: `save_img` image dimension is $n_dims - not saving")
    end
end

"""
    save_img(img::Array{UInt8,2}, filename::String)

Save an image to the input path.
"""
function save_img(img::Array{UInt8}, filename::String)
    mini, maxi = minimum(img), maximum(img)
    img = Float64.(img)
    if mini != maxi
        img = adjust_histogram(img, LinearStretching())
    end
    save_img(img, filename)
end

"""
    split_rgb(img::Array{UInt8,3})

Split input array in three channels.
"""
function split_rgb(img::Array{UInt8,3})
    @assert size(img)[1] == 3
    [img[i, :, :] for i in 1:3]
end


text:
import CPUSummary
# Increasing the number of threads must be done before importing Octavian
if Threads.nthreads() > 1 && Sys.CPU_THREADS > 1 && CPUSummary.num_cores() == 1
    CPUSummary.num_cores() = CPUSummary.static(2)
end

import Octavian

import Aqua
import BenchmarkTools
import ForwardDiff
import InteractiveUtils
import LinearAlgebra
import LoopVectorization
import Random
import Test
import VectorizationBase

using Test: @testset, @test, @test_throws, @inferred

include("test_suite_preamble.jl")

@info("VectorizationBase.num_cores() is $(VectorizationBase.num_cores())")
@info("Octavian.OCTAVIAN_NUM_TASKS[] is $(Octavian.OCTAVIAN_NUM_TASKS[]) tasks")

Random.seed!(123)
Octavian.debug() = true

include("block_sizes.jl")
include("init.jl")
include("integer_division.jl")
include("macrokernels.jl")
include("_matmul.jl")
coverage || include("matmul_main.jl")
include("matmul_coverage.jl")
include("utils.jl")
if sizeof(Int) >= 8 || !Sys.iswindows()
  include("forward_diff.jl")
end

include("aqua.jl") # run the Aqua.jl tests last


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
ENV["JULIA_DEBUG"] = SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

using Test

using Random
randseed = rand(UInt64)
@show randseed
Random.seed!(randseed)
#Random.seed!(0xdffe22467edd2cf0)  # needs more lenient SeparableArrays test

using LinearAlgebra

long_tests = true

## runtests:
@testset "SphericalFourierBesselDecompositions" begin
    include("test_separablearrays.jl")
    include("test_gnl.jl")
    include("test_cat2anlm.jl")
    include("test_windows.jl")
    include("test_theory.jl")  # should probably rename from "theory" to "lnn<->nlm conversion"
    include("test_modes.jl")
    include("test_wigner_chains.jl")
    include("test_nditerators.jl")
    include("test_window_chains.jl")
    include("test_covariance.jl")
    include("test_sfb.jl")
end

# vim: set sw=4 et sts=4 :


text:
using DataFrames
using Dates
using Plots
using InfrastructureSystemsPlots
using InfrastructureSystems
using Test
using TestSetExtensions
using Weave

const ISP = InfrastructureSystemsPlots
const IS = InfrastructureSystems
path = joinpath(pwd(), "plots")
!isdir(path) && mkdir(path)

function test_plots(file_path::String)
    include("test_data.jl")

    @testset "testing results sorting" begin
        Variables = Dict(:P_ThermalStandard => [:one, :two])
        sorted = ISP.sort_data(res; Variables = Variables)
        sorted_two = ISP.sort_data(res)
        sorted_names = [:one, :two]
        sorted_names_two = [:one, :three, :two]
        @test names(sorted.variables[:P_ThermalStandard]) == sorted_names
        @test names(sorted_two.variables[:P_ThermalStandard]) == sorted_names_two
    end

    @testset "testing bar plot" begin
        results =
            ISP.Results(res.variables, res.total_cost, res.optimizer_log, res.time_stamp)
        for name in keys(results.variables)
            variable_bar = ISP.get_bar_plot_data(results, string(name))
            sort = sort!(names(results.variables[name]))
            sorted_results = res.variables[name][:, sort]
            for i in 1:length(sort)
                @test isapprox(
                    variable_bar.bar_data[i],
                    sum(sorted_results[:, i]),
                    atol = 1.0e-4,
                )
            end
            @test typeof(variable_bar) == ISP.BarPlot
        end
        bar_gen = ISP.get_bar_gen_data(results)
        @test typeof(bar_gen) == ISP.BarGeneration
    end

    @testset "testing size of stack plot data" begin
        results =
            ISP.Results(res.variables, res.total_cost, res.optimizer_log, res.time_stamp)
        for name in keys(results.variables)
            variable_stack = ISP.get_stacked_plot_data(results, string(name))
            data = variable_stack.data_matrix
            legend = variable_stack.labels
            @test size(data) == size(res.variables[name])
            @test length(legend) == size(data, 2)
            @test typeof(variable_stack) == ISP.StackedArea
        end
        gen_stack = ISP.get_stacked_generation_data(results)
        @test typeof(gen_stack) == ISP.StackedGeneration
    end

    @testset "testing plot production" begin
        ISP.bar_plot(res; save = file_path, display = false)
        ISP.stack_plot(res; save = file_path, display = false)
        ISP.fuel_plot(res, generators; save = file_path, display = false)
        list = readdir(file_path)
        @test list == [
            "Bar_Generation.png",
            "Fuel_Bar.png",
            "Fuel_Stack.png",
            "P_RenewableDispatch_Bar.png",
            "P_RenewableDispatch_Stack.png",
            "P_ThermalStandard_Bar.png",
            "P_ThermalStandard_Stack.png",
            "Stack_Generation.png",
        ]
    end

    @testset "testing multi-plot production" begin
        ISP.bar_plot([res; res]; save = file_path, display = false)
        ISP.stack_plot([res; res]; save = file_path, display = false)
        ISP.fuel_plot([res; res], generators; save = file_path, display = false)
        list = readdir(file_path)
        @test list == [
            "Bar_Generation.png",
            "Fuel_Bar.png",
            "Fuel_Stack.png",
            "P_RenewableDispatch_Bar.png",
            "P_RenewableDispatch_Stack.png",
            "P_ThermalStandard_Bar.png",
            "P_ThermalStandard_Stack.png",
            "Stack_Generation.png",
        ]
    end

    @testset "testing report production" begin
        ISP.report(res, file_path)
        @test isfile(joinpath(file_path, "report_design.pdf"))
    end
end
try
    test_plots(path)
finally
    @info("removing test files")
    rm(path, recursive = true)
end


text:
""" Finite presentations of a model of a generalized algebraic theory (GAT).

We support two methods for defining models of a GAT: as Julia objects using the
`@instance` macro and as syntactic objects using the `@present` macro.
Instances are useful for casting generic data structures, such as matrices,
abstract tensor systems, and wiring diagrams, in categorical language.
Presentations define small categories by generators and relations and are useful
in applications like knowledge representation.
"""
module Present
export @present, Presentation, Equation, generator, generators, has_generator,
  equations, add_generator!, add_generators!, add_definition!, add_equation!,
  merge_presentation!

using Base.Meta: ParseError
using Compat
import DataStructures: OrderedSet
using Match

using ..Meta, ..Syntax
import ..Syntax: parse_json_sexpr, to_json_sexpr

# Data types
############

const GenExpr = GATExpr{:generator}
const Equation = Pair{<:GATExpr}{<:GATExpr}

mutable struct Presentation{T}
  generators::OrderedSet{GenExpr}
  generators_by_name::Dict{T,GenExpr}
  equations::OrderedSet{Equation}
end
Presentation(T::Type) = Presentation{T}(
  OrderedSet{GenExpr}(), Dict{T,GenExpr}(), OrderedSet{Equation}())
Presentation() = Presentation(Symbol)

# Presentation
##############

""" Get all generators of a presentation.
"""
function generators(pres::Presentation)::Vector
  collect(pres.generators)
end
function generators(pres::Presentation, typ::Type)::Vector
  filter(x -> isa(x,typ), collect(pres.generators))
end

""" Retrieve generators by name.
"""
function generator(pres::Presentation{T}, name) where T
  pres.generators_by_name[convert(T, name)]
end
function generators(pres::Presentation, names)::Vector
  [ generator(pres, name) for name in names ]
end

""" Does the presentation contain a generator with the given name?
"""
function has_generator(pres::Presentation{T}, name) where T
  haskey(pres.generators_by_name, convert(T, name))
end

""" Add a generator to a presentation.
"""
function add_generator!(pres::Presentation{T}, expr::GenExpr) where T
  name = first(expr)
  if name != nothing
    name = convert(T, name)
    if haskey(pres.generators_by_name, name)
      error("Name $name already defined in presentation")
    end
    pres.generators_by_name[name] = expr
  end
  push!(pres.generators, expr)
  return expr
end

""" Add multiple generators to a presentation.
"""
function add_generators!(pres::Presentation, exprs)
  for expr in exprs
    add_generator!(pres, expr)
  end
end

""" Get all equations of a presentation.
"""
function equations(pres::Presentation)::Vector
  collect(pres.equations)
end

""" Add an equation between terms to a presentation.
"""
function add_equation!(pres::Presentation, lhs::GATExpr, rhs::GATExpr)
  push!(pres.equations, lhs => rhs)
end

""" Add a generator defined by an equation.
"""
function add_definition!(pres::Presentation, name::Symbol, rhs::GATExpr)
  generator = Syntax.generator_like(rhs, name)
  add_generator!(pres, generator)
  add_equation!(pres, generator, rhs)
  return generator
end

""" Merge the second presentation into the first.

The first presentation is mutated and returned; the second is not.
"""
function merge_presentation!(pres::Presentation{T}, other::Presentation{T}) where T
  union!(pres.generators, other.generators)
  merge!(pres.generators_by_name, other.generators_by_name)
  union!(pres.equations, other.equations)
  return pres
end

# Presentation macro
####################

""" Define a presentation using a convenient syntax.
"""
macro present(head, body)
  name, syntax_name = @match head begin
    Expr(:call, [name::Symbol, arg::Symbol]) => (name, arg)
    _ => throw(ParseError("Ill-formed presentation header $head"))
  end
  let_expr = Expr(:let, Expr(:block), translate_presentation(syntax_name, body))
  esc(Expr(:(=), name, let_expr))
end

""" Translate a presentation in the DSL to a block of Julia code.
"""
function translate_presentation(syntax_name::Symbol, body::Expr)::Expr
  @assert body.head == :block
  code = Expr(:block)
  append_expr!(code, :(_presentation = $(module_ref(:Presentation))()))
  for expr in strip_lines(body).args
    append_expr!(code, translate_expr(syntax_name, expr))
  end
  append_expr!(code, :(_presentation))
  return code
end

""" Translate a single statement in the presentation DSL to Julia code.
"""
function translate_expr(syntax_name::Symbol, expr::Expr)::Expr
  @match expr begin
    Expr(:(::), [name::Symbol, type_expr]) =>
      translate_generator(syntax_name, name, type_expr)
    Expr(:(::), [type_expr]) =>
      translate_generator(syntax_name, nothing, type_expr)
    Expr(:(:=), [name::Symbol, def_expr]) =>
      translate_definition(name, def_expr)
    Expr(:(=), _) => expr
    Expr(:call, [:(==), lhs, rhs]) => translate_equation(lhs, rhs)
    _ => throw(ParseError("Ill-formed presentation statement $expr"))
  end
end

""" Translate declaration of a generator.
"""
function translate_generator(syntax_name::Symbol, name::Union{Symbol,Nothing},
                             type_expr)::Expr
  function rewrite(expr)
    @match expr begin
      Expr(:call, [name::Symbol, args...]) =>
        Expr(:call, GlobalRef(Syntax, :invoke_term), syntax_name,
             QuoteNode(name), map(rewrite, args)...)
      _ => expr
    end
  end
  type_name, args = @match type_expr begin
    name::Symbol => (name, [])
    Expr(:call, [name::Symbol, args...]) => (name, args)
    _ => throw(ParseError("Ill-formed type expression $type_expr"))
  end
  call_expr = Expr(
    :call, module_ref(:add_generator!), :_presentation,
    rewrite(Expr(:call, type_name,
                 isnothing(name) ? :nothing : QuoteNode(name), args...)))
  isnothing(name) ? call_expr : :($(name) = $call_expr)
end

""" Translate definition of generator in terms of other generators.
"""
function translate_definition(name::Symbol, def_expr)::Expr
  call_expr = Expr(:call, module_ref(:add_definition!), :_presentation,
                   QuoteNode(name), def_expr)
  :($name = $call_expr)
end

""" Translate declaration of equality between terms.
"""
function translate_equation(lhs, rhs)::Expr
  Expr(:call, module_ref(:add_equation!), :_presentation, lhs, rhs)
end

module_ref(sym::Symbol) = GlobalRef(Present, sym)

# Serialization
###############

function to_json_sexpr(pres::Presentation, expr::GATExpr)
  to_json_sexpr(expr;
    by_reference = name -> has_generator(pres, name))
end

function parse_json_sexpr(pres::Presentation, syntax_module::Module, sexpr)
  parse_json_sexpr(syntax_module, sexpr;
    parse_reference = name -> generator(pres, name))
end
function parse_json_sexpr(pres::Presentation{Symbol}, syntax_module::Module, sexpr)
  parse_json_sexpr(syntax_module, sexpr;
    symbols = true,
    parse_reference = name -> generator(pres, name))
end

end


text:
# Autogenerated wrapper script for ASL_jll for i686-linux-musl
export libasl

JLLWrappers.@generate_wrapper_header("ASL")
JLLWrappers.@declare_library_product(libasl, "libasl.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libasl,
        "lib/libasl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    supnorm(u, v)

Find the largest element-wise difference between two collections.
"""
function supnorm(u, v)
    return maximum(abs.(u .- v))
end

"""
    fixedpoint(iterate!, init[, cache]; distance=supnorm, maxiter=1000, tol=1e-10)

Repeatedly call iterate!, starting with `init`, until an approximate fixed point is reached.

Two iterations must be stored in order to calculate the distance between operations.
These iterations are stored in `cache`, which can be optionally supplied to avoid allocations.

# Examples
```jldoctest
julia> init = [1.0, 2.0];

julia> cache = (similar(init), similar(init));

julia> function iterate!(x1, x0, β)
           x1 .= β .* x0
       end;

julia> fixedpoint((x1, x0) -> iterate!(x1, x0, 0.9), init, cache)
(value = [4.165953662869907e-10, 8.331907325739814e-10], err = 9.257674806377564e-11, iter = 205, converged = true)
```
"""
function fixedpoint(iterate!, init, cache; distance=supnorm, maxiter=1000, tol=1e-10)
    x0, x1 = cache
    copy!(x0, init)
    err = zero(eltype(init))
    err_old = zero(eltype(init))

    for iter in 1:maxiter
        iterate!(x1, x0)
        err = distance(x1, x0)
        if iter > 1 && err >= err_old
            @warn "Iteration: $iter; Error increased by $((err / err_old - 1) * 100)%"
        end
        @debug "Iteration: $iter; Current error: $err"
        if err < tol
            return (value=x1, err=err, iter=iter, converged=true)
        else
            copy!(x0, x1)
            err_old = err
        end
    end
    return (value=x1, err=err, iter=maxiter, converged=false)
end

function fixedpoint(iterate!, init; distance=supnorm, maxiter=1000, tol=1e-10)
    cache = (copy(init), copy(init))
    return fixedpoint(iterate!, init, cache; distance=distance, maxiter=maxiter, tol=tol)
end


text:
#populate array
# %%

# Atom Editorで起動 uber-junoがすでに導入されているとする。
# ここら辺にカーソルを移動
# option + enter を押す。
# 次の`# %%` が出てくるまでの部分まで実行される
println("ここから処理が始まるよ")
arr = [1,2,3]
@show arr

arr = Float64[]
push!(arr,1)
@show arr

arr = Int64[]
push!(arr,1)
@show arr

println("ここで処理が終わるよ")

# %%
@show collect(1:7)
@show size(collect(1:7),1)

arr = Float64[]
sizehint!(arr,6)
@show size(arr)
for i in 1:10
    push!(arr,i)
end

@show arr

# %%

arr = collect(1:10)
@show arr[3:7]

# %%
@show println(Array{Any}(undef,4))
@show fill!(arr,42)

# %% concatenate arrays

a₁ = [-1,-2,-3]
a₂ = [4,5,6]

append!(a₁,a₂)
@show a₁

# you can't do
#append!(a₁,9)

# %%

arr = collect(1:8)
@show arr
p=pop!(arr)
@show p
@show arr

p=popfirst!(arr)
@show p
@show arr
pushfirst!(arr,-1)
@show arr

# %%
# splice

arr = collect(1:8)

@show arr
splice!(arr,3)
@show arr

# %%

#=
Checking wheather an array contains an element
=#
arr = collect(1:8)

@show　4 in arr
@show -1 in arr
@show in(4,arr)
@show in(-1,arr)

#=
julia> import LinearAlgebra:dot

julia> f = x -> dot(x,x); g = x -> sum(x .* x);

julia> arr =collect(1:1000);

julia> @benchmark f(arr)
BenchmarkTools.Trial:
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     264.665 ns (0.00% GC)
  median time:      268.938 ns (0.00% GC)
  mean time:        326.445 ns (10.09% GC)
  maximum time:     276.418 μs (99.82% GC)
  --------------
  samples:          10000
  evals/sample:     328

julia> @benchmark g(arr)
BenchmarkTools.Trial:
  memory estimate:  7.95 KiB
  allocs estimate:  2
  --------------
  minimum time:     715.370 ns (0.00% GC)
  median time:      1.134 μs (0.00% GC)
  mean time:        1.772 μs (37.87% GC)
  maximum time:     897.657 μs (99.85% GC)
  --------------
  samples:          10000
  evals/sample:     108

julia> f(arr) == g(arr)

=#

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Plots, NBodyIPs, JuLIP
using NBodyIPs: Cutoff
using BenchmarkTools

##
r0, rnn, rc = 0.0, 1.0, 2+√2
xx = range(r0-0.5, rc+0.5, length=300)

ξ1 = @analytic r -> (r - rnn)/(rc-r0) * ((r-r0)/(rc-rnn) + (r-rc)/(r0-rnn))
ξ2 = @analytic r -> (r-rnn) * exp(-r)

plot(xx, ξ1.(xx))
scatter!([r0, rnn, rc], [-1,0,1], ms=10)

##

using LinearAlgebra, Roots
r0, rnn, rc = 0.75, 1.0, 2.1
xx = range(r0-0.5, rc+0.5, length=300)

# ξ(r) = A (1/(1+r) - 1/2) + B (r-1)
c = qr( [1/(r0+1)-1/2  r0-1; 1/(rc+1)-1/2  rc-1] ) \ [-1.0, 1.0]
ξ2 = r -> c[1]*(1/(r+1)-1/2) + c[2] * (r - 1)

# ξ(r) = C(exp( λ (r/r0 - 1) ) - 1)
gg = λ -> exp( λ*(rc/rnn - 1) ) + exp( λ*(r0/rnn - 1) ) - 2
@show λopt = find_zero(gg, -2.5)
C = 1 / (exp( λopt*(rc/rnn - 1) ) - 1)
ξ3  = r -> C * (exp( λopt*(r/rnn - 1) ) - 1)

# cubic, minimising curvature
#  (r-rnn)/(rc-rnn) * (1 + A (r-rc) + B (r-rc)^2)
B = (7*(r0^2 + r0*rc + rc^2 - 3*(r0 + rc)*rnn + 3*rnn^2)) / (
     2*(r0 - rc)*(r0 - rnn)*(2*r0^2 + 3*r0*rc + 2*rc^2 -
     7*(r0 + rc)*rnn + 7*rnn^2))
A = 2/(-r0 + rc) + 3/(4*(-r0 + rnn)) + (7*(rc - rnn)) / (
    4*(2*r0^2 + 3*r0*rc + 2*rc^2 - 7*(r0 + rc)*rnn + 7*rnn^2))
ξ4 = r -> (r-rnn)/(rc-rnn) * (1 + A*(r-rc) + B*(r-rc)^2)
χ = x -> (x+1)^2 * (x-1)^2

# another polynomial attempt...

plot(xx, ξ1.(xx))
plot!(xx, ξ2.(xx))
plot!(xx, ξ3.(xx))
plot!(xx, ξ4.(xx))
scatter!([r0, rnn, rc], [-1,0,1], ms=10)

##
χ = x -> (x+1)^2 * (x-1)^2
penv2s = r -> χ(ξ3(r)) * (r0 < r < rc)
plot(xx, penv2s)
scatter!([r0, rc], [0,0])

##
cutcos = Cutoff((:cos2s, 0.7, 0.9, 2.3, 3.0))
cutpenv = Cutoff((:penv2s, 2, 0.7, 1.0, 3.0))

r0, rnn, rc = 0.7, 1.0, 3.0
ξ = @analytic r -> (r - rnn)/(rc-r0) * ((r-r0)/(rc-rnn) + (r-rc)/(r0-rnn))

env = @analytic x -> (x-1)^2 * (x+1)^2
f = r -> env(ξ(r)) * (r0 < r < rc)

xx = range(0, 3.5, length=300)
P1 = plot(xx, cutcos.f.(xx), label = "cos2s")
plot!(xx, cutpenv.f.(xx), label = "penv2s")

## one-sided cutoffs
cosenv(x) = (1 + cos(pi*x))/2
p2env(x) = (1-x)^2
p2envsym(x) = (1-x)^2*(1+x)^2
p3env(x) = (1-x)^3
p3envsym(x) = (1-x)^3*(1+x)^3

xx = range(0, 1, length=100)
P1 = plot(xx, cosenv.(xx), label = "cosenv")
plot!(xx, p2env.(xx), label = "p2env")
plot!(xx, p2envsym.(xx), label = "p2envsym")
plot!(xx, p3env.(xx), label = "p3env")
plot!(xx, p3envsym.(xx), label = "p3envsym")

## two-sided cutoffs

xx = range(-1, 1, length=100)
P2 = plot(xx, cosenv.(xx), label = "cosenv")
plot!(xx, p2envsym.(xx), label = "p2envsym")
plot!(xx, p3envsym.(xx), label = "p3envsym")
plot!(xx, cosenv.(xx).^2, label = "cosenv^2")

plot(P1, P2)

##
using Roots

struct PEnv2s{T}
   C::T
   λ::T
   rnn::T
end

struct PEnv2sB{T}
   C::T
   λ::T
   rnn::T
end

function PEnv2s(r0, rnn, rc)
   gg = λ -> exp( λ*(rc/rnn - 1) ) + exp( λ*(r0/rnn - 1) ) - 2
   @show λopt = find_zero(gg, -2.5)
   C = 1 / (exp( λopt*(rc/rnn - 1) ) - 1)
   return PEnv2s(C, λopt)
end

Base.@pure function (p::PEnv2s)(r)
   rnn = 1.0
   x = @fastmath(p.C * (exp( p.λ*(r - 1) ) - 1))
   a = x*x-1
   return a*a
end

x = 1.456

function runN(f, x, N)
   s = 0.0
   for n = 1:N
      x += 1e-9
      s += f(x)
   end
   return s
end

penv = Cutoff((:penv2s, 2, 0.7, 1.0, 2.5))
penv2 = PEnv2s(0.7, 1.0, 2.5)
cos2s = Cutoff((:cos2s, 0.7, 0.9, 1.7, 2.5))

@btime(runN($(penv.f), $x, 1_000_000))
@btime(runN($(cos2s.f), $x, 1_000_000))
@btime(runN($(penv2), $x, 1_000_000))

@btime runN($cosenv, $x, 1_000)
@btime runN($p2env, $x, 1_000)
@btime runN($p2envsym, $x, 1_000)
@btime runN($p3env, $x, 1_000)
@btime runN($p3envsym, $x, 1_000)

macro horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :( $(esc(p[i])) + t * $ex )
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

f1(r) = 2.345 * (exp( 2.345*r - 1.234) - 1.213)
f2(r) = (@horner(r, 1.234, 3.234, 4.234, 4.234, 4.234, 3.2314, 4.234))
@btime runN($f1, $x, 1_000)
@btime runN($f2, $x, 1_000)
@btime runN($exp, $x, 1_000)
@btime runN($cos, $x, 1_000)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# test implementation
struct NewTrait<:NewType
    name::Symbol
    super::Any
end
Base.show(io::IO, x::NewTrait) = print(io, x.name)
TypeEmulator.supertypes(x::NewTrait) = x.super
TypeEmulator.is_pushed(::Type{NewTrait}) = true
isnewsub(a::NewTrait, b::NewTrait) = a != b && b ∈ linearize(a)

# Test

Anny = NewTrait(:Any, NewTrait[])

Animal = NewTrait(:Animal, [Anny])
HasLegs = NewTrait(:HasLegs, [Animal])
FourLegged = NewTrait(:FourLegged, [HasLegs])
TwoLegged = NewTrait(:TwoLegged, [HasLegs])
Furry = NewTrait(:Furry, [Animal])
Cat = NewTrait(:Cat, [Animal, FourLegged, Furry])

@test linearize(Anny) == [Anny]
@test linearize(Animal) == reverse([Animal, Anny])
@test linearize(FourLegged) == reverse([FourLegged, HasLegs, Animal, Anny])
@test linearize(Cat) == reverse([Cat, Furry, FourLegged, HasLegs, Animal, Anny])

@test isnewsub(Cat, Animal)
@test isnewsub(Cat, HasLegs)
@test !isnewsub(Anny, Anny)
@test !isnewsub(Cat, Cat)
@test isnewsub(FourLegged, HasLegs)


text:
function scores(a::AbstractArray; latest = false, top=0)
    if top > 0
        top > length(a) ? (return reverse(sort(a))) : nothing
        if top == 1
            return sort(a)[end]
        end
        return reverse(sort(a)[end-top+1:end])
    end

    latest ? (return a[end]) : return a
end

max_stars_repo_name: josobar/julia

max_issues_repo_name: josobar/julia

max_forks_repo_name: josobar/julia


text:
function pmttest()
	#@testset "pmt_test.jl" begin
        fword  = "/Users/jj/JuliaProjects/LaserLab/test/C1--Trace--00010.trc"
        fbyte  = "/Users/jj/JuliaProjects/LaserLab/test/C1--Trace--00108.trc"

        iob = open(fbyte, "r")
        iow = open(fword, "r")

        @test typeof(iob) == IOStream
        @test typeof(iow) == IOStream

        WAVEDESCb = LaserLab.wavedesc(iob)
        WAVEDESCw = LaserLab.wavedesc(iow)
        @test WAVEDESCb == WAVEDESCw

        @test LaserLab.readword(iob, WAVEDESCb + 32) == 0  #Byte
        @test LaserLab.readword(iow, WAVEDESCw + 32) == 1  #Word

        @test LaserLab.readlong(iow, WAVEDESCw + 60) == 2 * LaserLab.readlong(iob, WAVEDESCb + 60)

        @test LaserLab.readstring(iow, WAVEDESCw + 76, 14) == "LECROYWS4104HD"
        @test LaserLab.readstring(iob, WAVEDESCb + 76, 14) == "LECROYWS4104HD"

        @test LaserLab.readfloat(iow, WAVEDESCw + 160) == LaserLab.readfloat(iob, WAVEDESCb + 160)

        @test LaserLab.readdouble(iow, WAVEDESCw+ 180) ≈ LaserLab.readdouble(iob, WAVEDESCb + 180)

        tsw = LaserLab.readtimestamp(iow, WAVEDESCw + 296)
        @test tsw.year == 2022
    #end
end

max_stars_repo_name: jjgomezcadenas/LaserLab

max_issues_repo_name: jjgomezcadenas/LaserLab

max_forks_repo_name: jjgomezcadenas/LaserLab


text:
# low-level wrappers

# grid graph
GCoptimizationGridGraph(width::GCOSiteID, height::GCOSiteID, numLabels::GCOLabelID) = @cxxnew GCoptimizationGridGraph(width, height, numLabels)
setSmoothCostVH(gco::Cxx.CppPtr, smoothArrayPtr::Ptr{GCOEnergyTermType}, vCostsPtr::Ptr{GCOEnergyTermType}, hCostsPtr::Ptr{GCOEnergyTermType})::Void = @cxx gco -> setSmoothCostVH(smoothArrayPtr, vCostsPtr, hCostsPtr)

# general graph
# TODO: support setAllNeighbors
GCoptimizationGeneralGraph(num_sites::GCOSiteID, num_labels::GCOLabelID) = @cxxnew GCoptimizationGeneralGraph(num_sites, num_labels)
setNeighbors(gco::Cxx.CppPtr, site1::GCOSiteID, site2::GCOSiteID, weight::GCOEnergyTermType=GCOEnergyTermType(1))::Void = @cxx gco -> setNeighbors(site1, site2, weight)

# expansion
# "Peforms expansion algorithm. Runs the number of iterations specified by max_num_iterations.
#  If no input specified, runs until convergence. Returns total energy of labeling."
expansion(gco::Cxx.CppPtr, max_num_iterations::Cint=Cint(-1))::GCOEnergyType = @cxx gco -> expansion(max_num_iterations)

# "Peforms  expansion on one label, specified by the input parameter alpha_label."
alpha_expansion(gco::Cxx.CppPtr, alpha_label::GCOLabelID)::GCObool = @cxx gco -> alpha_expansion(alpha_label)

# swap
# "Peforms swap algorithm. Runs it the specified number of iterations.
#  If no input is specified,runs until convergence."
swap(gco::Cxx.CppPtr, max_num_iterations::Cint=Cint(-1))::GCOEnergyType = @cxx gco -> swap(max_num_iterations)

# "Peforms  swap on a pair of labels, specified by the input parameters alpha_label, beta_label"
alpha_beta_swap(gco::Cxx.CppPtr, alpha_label::GCOLabelID, beta_label::GCOLabelID)::Void = @cxx gco -> alpha_beta_swap(alpha_label, beta_label)

# "Peforms swap on a pair of labels, specified by the input parameters alpha_label, beta_label,
#  only on the sitess in the specified arrays, alphaSites and betaSitess, and the array sizes
#  are, respectively, alpha_size and beta_size"
alpha_beta_swap(gco::Cxx.CppPtr, alpha_label::GCOLabelID, beta_label::GCOLabelID, alphaSitesPtr::Ptr{GCOSiteID}, alpha_size::GCOSiteID, betaSitesPtr::Ptr{GCOSiteID}, beta_size::GCOSiteID)::Void = @cxx gco -> alpha_beta_swap(alpha_label, beta_label, alphaSitesPtr, alpha_size, betaSitesPtr, beta_size)

# setDataCost
# "Set cost for all (SiteID,LabelID) pairs. Default data cost is all zeros."
function setDataCost(gco::Cxx.CppPtr, dataCostFunc::Function)::Void
    const fn = cfunction(dataCostFunc, GCOEnergyTermType, (GCOSiteID, GCOLabelID))
    @cxx gco -> setDataCost(icxx"($GCOEnergyTermType (*)($GCOSiteID,$GCOLabelID))($fn);")
end
function setDataCost(gco::Cxx.CppPtr,  dataCostFuncExtra::Function, extraData::Ptr{Void})::Void
    const fn = cfunction(dataCostFuncExtra, GCOEnergyTermType, (GCOSiteID, GCOLabelID, Ptr{Void}))
    @cxx gco -> setDataCost(icxx"($GCOEnergyTermType (*)($GCOSiteID,$GCOLabelID,void *))($fn);", extraData)
end
setDataCost(gco::Cxx.CppPtr, dataArrayPtr::Ptr{GCOEnergyTermType})::Void = @cxx gco -> setDataCost(dataArrayPtr)
setDataCost(gco::Cxx.CppPtr, s::GCOSiteID, l::GCOLabelID, e::GCOEnergyTermType)::Void = @cxx gco -> setDataCost(s, l, e)

# TODO: support setDataCostFunctor

# "Set cost of assigning 'l' to a specific subset of sites.
#  The sites are listed as (SiteID,cost) pairs."
setDataCost(gco::Cxx.CppPtr, l::GCOLabelID, costsPtr::Ptr{GCOSparseDataCost}, count::GCOSiteID)::Void = @cxx gco -> setDataCost(l, costsPtr, count)

# setSmoothCost
# "Set cost for all (LabelID,LabelID) pairs; the actual smooth cost is then weighted
#  at each pair of on neighbors. Defaults to Potts model (0 if l1==l2, 1 otherwise)"
function setSmoothCost(gco::Cxx.CppPtr, smoothCostFunc::Function)::Void
    const fn = cfunction(smoothCostFunc, GCOEnergyTermType, (GCOSiteID, GCOSiteID, GCOLabelID, GCOLabelID))
    @cxx gco -> setSmoothCost(icxx"($GCOEnergyTermType (*)($GCOSiteID,$GCOSiteID,$GCOLabelID,$GCOLabelID))($fn);")
end
function setSmoothCost(gco::Cxx.CppPtr, smoothCostFuncExtra::Function, extraData::Ptr{Void})::Void
    const fn = cfunction(smoothCostFuncExtra, GCOEnergyTermType, (GCOSiteID, GCOSiteID, GCOLabelID, GCOLabelID, Ptr{Void}))
    @cxx gco -> setSmoothCost(icxx"($GCOEnergyTermType (*)($GCOSiteID,$GCOSiteID,$GCOLabelID,$GCOLabelID,void *))($fn));", extraData)
end
setSmoothCost(gco::Cxx.CppPtr, l1::GCOLabelID, l2::GCOLabelID, e::GCOEnergyTermType)::Void = @cxx gco -> setSmoothCost(l1, l2, e)
setSmoothCost(gco::Cxx.CppPtr, smoothArrayPtr::Ptr{GCOEnergyTermType})::Void = @cxx gco -> setSmoothCost(smoothArrayPtr)

# TODO: support setSmoothCostFunctor

# setLabelCost
# "Sets the cost of using label in the solution.
#  Set either as uniform cost, or an individual per-label cost."
setLabelCost(gco::Cxx.CppPtr, cost::GCOEnergyTermType)::Void = @cxx gco -> setLabelCost(cost)
setLabelCost(gco::Cxx.CppPtr, costArrayPtr::Ptr{GCOEnergyTermType})::Void = @cxx gco -> setLabelCost(costArrayPtr)
setLabelSubsetCost(gco::Cxx.CppPtr, labelsPtr::Ptr{GCOLabelID}, numLabels::GCOLabelID, cost::GCOEnergyTermType)::Void = @cxx gco -> setLabelSubsetCost(labelsPtr, numLabels, cost)

# whatLabel
# "Returns current label assigned to input site"
whatLabel(gco::Cxx.CppPtr, site::GCOSiteID)::GCOLabelID = @cxx gco -> whatLabel(site)
whatLabel(gco::Cxx.CppPtr, start::GCOSiteID, count::GCOSiteID, labelingPtr::Ptr{GCOLabelID})::Void = @cxx gco -> whatLabel(start, count, labelingPtr)

# setLabel
# "This function can be used to change the label of any site at any time"
setLabel(gco::Cxx.CppPtr, site::GCOSiteID, label::GCOLabelID)::Void = @cxx gco -> setLabel(site, label)

# setLabelOrder
# "setLabelOrder(false) sets the order to be not random; setLabelOrder(true)
#  sets the order to random. By default, the labels are visited in non-random order
#  for both the swap and alpha-expansion moves
#  Note that srand() must be initialized with an appropriate seed in order for
#  random order to take effect!"
setLabelOrder(gco::Cxx.CppPtr, isRandom::GCObool)::Void = @cxx gco -> setLabelOrder(isRandom)
setLabelOrder(gco::Cxx.CppPtr, orderPtr::Ptr{GCOLabelID}, size::GCOLabelID)::Void = @cxx gco -> setLabelOrder(orderPtr, size)

# compute_energy
# "Returns total energy for the current labeling"
compute_energy(gco::Cxx.CppPtr)::GCOEnergyType = @cxx gco -> compute_energy()

# give*energy
# "Returns separate Data, Smooth, and Label energy of current labeling"
giveDataEnergy(gco::Cxx.CppPtr)::GCOEnergyType = @cxx gco -> giveDataEnergy()
giveSmoothEnergy(gco::Cxx.CppPtr)::GCOEnergyType = @cxx gco -> giveSmoothEnergy()
giveLabelEnergy(gco::Cxx.CppPtr)::GCOEnergyType = @cxx gco -> giveLabelEnergy()

# numSites & numLabels
# "Returns number of sites/labels as specified in the constructor"
numSites(gco::Cxx.CppPtr)::GCOSiteID = @cxx gco -> numSites()
numLabels(gco::Cxx.CppPtr)::GCOLabelID = @cxx gco -> numLabels()

# setVerbosity
# "Prints output to stdout during exansion/swap execution.
#   0 => no output
#   1 => cycle-level output (cycle number, current energy)
#   2 => expansion-/swap-level output (label(s), current energy)"
setVerbosity(gco::Cxx.CppPtr, level::Cint)::Void = @cxx gco -> setVerbosity(level)


text:
mutable struct ZTMetaPopulation{T<:LanguageDynamics.AbstractInteractor}
  # Census (of local populations)
  census::Array{LanguageDynamics.PoolPopulation{T}}

  # Lattice size (side length)
  size::Int

  # Distance decay exponent
  de::Float64

  # Probability mass function for movement
  probabilities::Array{Float64}

  # Distance function
  distance::Function

  # Range of possible distances; will be different for different cells since
  # we assume fixed lattice boundaries
  possible_distances::Dict{Int, Array{Int}}

  # How many actions have occurred
  nactions::Int

  # Cells at a given distance from given cell
  lookup::Dict{Tuple{Int, Int}, Array{Int}}
end

function ZTMetaPopulation{T}(K::Int, lambda::Float64) where {T<:LanguageDynamics.AbstractInteractor}
  # Possible distances on lattice
  possdist = collect(0:2*(K-1))

  # Possible distances for each cell; initialized here as empty
  possdistcell = Dict{Int, Array{Int}}()

  # Metapopulation; array of local populations which are pools
  mp = Array{LanguageDynamics.PoolPopulation{T}}(undef, 0)

  # Construct lookup table and add local populations. This is
  # intensive, but will only happen once (when the metapopulation
  # is constructed).
  lup = Dict{Tuple{Int, Int}, Array{Int}}()
  cellID1 = 0
  for x in 1:K
    for y in 1:K
      cellID1 += 1
      point1 = LanguageDynamics.Point{Int}(x,y)

      push!(mp, LanguageDynamics.PoolPopulation{T}())

      for dist in possdist
        cellID2 = 0
        for i in 1:K
          for j in 1:K
            cellID2 += 1
            point2 = LanguageDynamics.Point{Int}(i,j)

            if LanguageDynamics.manhattan(point1, point2) == dist
              # if this cell does not yet have an entry, initialize
              if !haskey(possdistcell, cellID1)
                possdistcell[cellID1] = Array{Int}(undef, 0)
              end

              # push this distance to the array
              push!(possdistcell[cellID1], dist)

              # if this pair does not yet exist in lookup, initialize
              if !haskey(lup, (cellID1, dist))
                lup[(cellID1, dist)] = Array{Tuple{Int, Int}}(undef, 0)
              end

              # push cell 2 to the lookup for this pair
              push!(lup[(cellID1, dist)], cellID2)
            end
          end
        end
      end

      # we only need cell IDs once in the array
      possdistcell[cellID1] = unique(possdistcell[cellID1])
    end
  end

  # Construct
  ZTMetaPopulation(mp, K, lambda, LanguageDynamics.zipf(2*(K-1) + 1, lambda), LanguageDynamics.manhattan, possdistcell, 0, lup)
end

function inject!(x::LanguageDynamics.AbstractInteractor, y::ZTMetaPopulation)
  # random local population
  rp = rand(1:length(y.census))

  # assign this to speaker's home, using speaker's meta variable
  x.meta = rp

  # inject into local population
  LanguageDynamics.inject!(x, y.census[rp])
end

function rendezvous!(x::ZTMetaPopulation)
  rend = LanguageDynamics.rendezvous!(rand(x.census))
  if rend != false
    x.nactions += 2
  end
end

function travel!(x::ZTMetaPopulation)
  # select random local population
  rp = rand(x.census)

  # select random speaker; fail if local population empty
  while length(rp.census) == 0
    rp = rand(x.census)
  end
  #length(rp.census) > 0 || return false
  rs_coor = rand(1:length(rp.census)) # NB: coordinate in *local* population
  rs = rp.census[rs_coor]

  # travel
  home = rs.meta
  possdist = x.possible_distances[home]
  dist = Distributions.wsample(possdist, LanguageDynamics.zipf(maximum(possdist) + 1, x.de))
  destination = rand(x.lookup[(home, dist)])
  LanguageDynamics.eject!(rs_coor, rp)
  LanguageDynamics.inject!(rs, x.census[destination])
end


text:
function vandermonde_matrix(offsets::Array{Int64, 1})
    num_points = size(offsets)[1]
    ret_matrix = reshape([0 // 1 for i in 1:num_points ^ 2],
                         (num_points, num_points))
    for i in 1:num_points
        ret_matrix[1, i] = 1 // 1
    end
    for i in 2:num_points
        for j in 1:num_points
            ret_matrix[i, j] = ret_matrix[i - 1, j] * offsets[j]
        end
    end
    return ret_matrix
end # vandermonde_matrix

function scheme_coefficients(offsets::Array{Int64, 1})
    num_points = size(offsets)[1]
    coefficient_matrix = vandermonde_matrix(offsets)
    right_vector = reshape([0 // 1 for i in 1:num_points],
                           (num_points, 1))
    right_vector[2, 1] = 1 // 1
    return reshape(coefficient_matrix \ right_vector, (num_points, ))
end # scheme_coefficients

function rational_flux_coefficients(offsets::Array{Int64, 1})
    num_points = size(offsets)[1]
    coeff = scheme_coefficients(offsets)
    ret = [0 // 1 for i in 1:num_points]
    for i in 2:num_points
        ret[i] = ret[i - 1] - coeff[i - 1]
    end
    return ret[2:num_points]
end # rational_flux_coefficients

function float_flux_coefficients(offsets::Array{Int64, 1})
    rational_coeff = rational_flux_coefficients(offsets)
    num_points = size(rational_coeff)[1]
    return [float(rational_coeff[i]) for i in 1:num_points]
end # float_flux_coefficients


text:
#
# This module is an entry point of jevo application. jevo is a
# Julia(http://julialang.org/) language based, digital organisms
# evolution simulator.It's used for study the evolutionary biology
# of self-replicating and evolving computer programs (digital
# organisms (https://en.wikipedia.org/wiki/Digital_organism)).
# This project similar to Avida(https://en.wikipedia.org/wiki/Avida),
# but works with more abstract language (Julia) instead of assembler.
# More details on blog(https://jevosite.wordpress.com/) and youtube
# (https://www.youtube.com/playlist?list=PL1NiKjXMaBimPuybPIXkVuO1MYy53XcdW)
# channel.
#
# jevo is also a client\server application. When you run one Manager
# (current file) we are talking about one Manager's instance. It's
# server and client at the same time, because other instances or remote
# terminals\visualizers applications may be connected to it. Every
# Manager manages virtual 2D world, runs virtual processor, parallel
# organisms codes and so on.
#
# You may run the manager in different ways. The simplies way is:
#   >julia src\AppSatellite.jl src\AppManager.jl
# More useful way will be the following:
#   >julia --color=yes src\AppSatellite.jl src\AppManager.jl modeQuiet
# It colors terminal output and reduces amount of default messages.
# Also, jevo supports many command line arguments like: recover,
# about, version and so on... You may use them at the end of command
# line. For example, to get current jevo version you have to type:
#   >julia src\AppManager.jl version
#
#
# TODO: describe "recover" command line argument
#
include("global/ImportFolders.jl")

import Config.@if_profile
import Config.@if_not_profile
import Config
import Manager
import Helper
import Backup
import CommandLine
import About
@if_profile import ProfileView
#
# Name of the command line argument, which tells the application
# to recover itself from last backup.
# TODO: command line commands should be moved to separate module
const ARG_RECOVER = "recover"
const ARG_ABOUT   = "about"
const ARG_INFO    = "info"
const ARG_VERSION = "version"
#
# This function starts the manager, world, organisms, server etc...
# It checks "recover" argument for recovering from last backup or
# runs in a common mode.
#
function main()
  local args::Dict{String, String} = CommandLine.create()
  local exitCode::Int
  local command::Pair{String, Function}
  local commands::Dict{String, Function} = Dict{String, Function}(
    ARG_ABOUT   => _onAbout,
    ARG_INFO    => _onAbout,
    ARG_VERSION => _onVersion,
    ARG_RECOVER => _onRecover
  )

  for command in commands
    if CommandLine.has(args, command[1])
      exitCode = command[2]()
      @if_profile return exitCode
      exit(exitCode)
    end
  end

  exitCode = _onRun()
  @if_profile return exitCode
  exit(exitCode)
end
#
# Just returns about info
#
function _onAbout()
  Helper.info(About.description())
  0
end
#
# Returns jevo version
#
function _onVersion()
  Helper.info(About.version())
  0
end
#
# Runs jevo application, but before recovers itself from
# last backup file in "backup" folder.
#
function _onRecover()
  local man::ManagerTypes.ManagerData
  local exitCode::Int
  local cfg::Config.ConfigData = Config.create(false)

  if cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Starting jevo...") end
  man = Manager.create(cfg)
  Manager.recover(man)
  #
  # Overrides recovered config by command line arguments
  #
  man.cfg.cmdLineArgs = CommandLine.create()
  man.cfg = Config.merge(man.cfg)
  exitCode = Int(!Manager.run(man, true)) # 1 - error, 0 - okay
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Quit jevo...") end

  exitCode
end
#
# Handles running by "default". Runs jevo application, starts
# Manager and evolution of digital organisms
#
function _onRun()
  local man::ManagerTypes.ManagerData
  local exitCode::Int
  local cfg::Config.ConfigData = Config.create(false)

  man = Manager.create(cfg)
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Starting jevo...") end
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Running from scratch...") end
  exitCode = Int(!Manager.run(Manager.create())) # 1 - error, 0 - okay
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Quit jevo...") end

  exitCode
end
#
# Application entry point. In case of profiling mode, application
# will be run in special profiling mode using ProfileView package.
# In this case fire chart will be drown at the end
#
@if_profile Profile.clear()
@if_profile @profile main()
@if_profile ProfileView.view()
@if_not_profile main()

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id="input", value = "initial value"),
    dcc_input(id="input2", value = "prevented"),
    html_div("initial",id="output")
end
callback!(app,
    Output("input2","value"),
    Input("input","value"),
    prevent_initial_call = true
    ) do value
    return value
end
callback!(app,
    Output("output","children"),
    Input("input2","value")
    ) do value
    return value
end

run_server(app)


text:
module UKE

using Reexport
@reexport using Kinetic
@reexport using Flux
@reexport using DiffEqFlux
@reexport using OrdinaryDiffEq
using Optim
using FileIO
using JLD2
using Plots
import Tracker

export device
export track, untrack, tracker_mode, zygote_mode
export vis_train
export ube_dfdt, ube_dfdt!
export step_ube!

include("widget.jl")
include("tracker.jl")
include("train.jl")
include("collision.jl")
include("solver.jl")

end # module


text:
# Check type stability of ODE_DEFAULT_NORM with mixed units and with ArrayPartition
using Test
using MechGlueDiffEqBase
using MechanicalUnits: @import_expand, dimension, NoDims, ∙
import MechanicalUnits: g, g⁻¹
import MechanicalUnits.Unitfu.numtype
@import_expand(km, N, s, m, km, kg, °, inch)
using DiffEqBase, OrdinaryDiffEq

@testset "Inferrable norm Unitless ArrayPartition" begin
    r0ul = [1131.340, -2282.343, 6672.423]
    v0ul = [-5.64305, 4.30333, 2.42879]
    rv0ul = @inferred ArrayPartition(r0ul,v0ul)
    @test @inferred(ODE_DEFAULT_NORM(rv0ul, 0.0)) === 2915.770473301504
end
@testset "Inferrable norm  ArrayPartition units" begin
    r0 = [1131.340, -2282.343, 6672.423]∙km
    v0 = [-5.64305, 4.30333, 2.42879]∙km/s
    rv0 = @inferred ArrayPartition(r0, v0)
    @test @inferred(ODE_DEFAULT_NORM(rv0, 0.0s)) === 41.02536038842316
end
@testset "Inferrable norm  ArrayPartition mixed units" begin
    r0 = [1.0km, 2.0km, 3.0m/s, 4.0m/s]
    v0 = [1.0km/s, 2.0km/s, 3.0m/s², 4m/s²]
    rv0 = @inferred ArrayPartition(r0, v0)
    @test @inferred(ODE_DEFAULT_NORM(rv0, 0.0s)) === 1.5811388300841898
end

@testset "Inferrable zero Unitless ArrayPartition" begin
    r0ul = [1131.340, -2282.343, 6672.423]
    v0ul = [-5.64305, 4.30333, 2.42879]
    rv0ul = ArrayPartition(r0ul, v0ul)
    @test @inferred(zero(rv0ul)) == ArrayPartition([0.0, 0.0, 0.0], [0.0, 0.0, 0.0])
    @test ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}} == typeof(zero(rv0ul))
end
@testset "Inferrable zero compatible units" begin
    r0 = [1131.340, -2282.343, 6672.423]∙km
    r1 = [1km, 2.0m]
    @test @inferred(zero(r0)) == [0.0, 0.0, 0.0]km
    @test @inferred(zero(r1)) == [0.0, 0.0]km
    rv0 = ArrayPartition(r0)
    @test @inferred(zero(rv0)) == [0.0, 0.0, 0.0]km
    @test typeof(zero(rv0)) == typeof(rv0)
    rv1 = ArrayPartition(r1)
    @test @inferred(zero(rv1)) == [0.0, 0.0]km
    @test typeof(zero(rv0)) == typeof(rv0)
end
@testset "Inferrable zero vector and ArrayPartition mixed units" begin
    r0 = [1.0km, -2km, 3m/s, 4m/s]
    @test @inferred(zero(r0)) == [0.0km, 0.0km, 0.0m/s, 0.0m/s]
    rv0 = @inferred ArrayPartition(r0)
    zer = @inferred zero(rv0)
    @test zer == [0.0km, 0.0km, 0.0m/s, 0.0m/s]
    @test typeof(zer) === typeof(rv0)
end

@testset "Inferrable similar vector and ArrayPartition compatible units" begin
    r0 = [1131.340, -2282.343, 6672.423]∙km
    simi = @inferred(similar(r0))
    @test all(typeof.(r0) == typeof.(simi))
    @test simi !== r0
    rv0 = @inferred ArrayPartition(r0)
    sima = @inferred similar(rv0)
    @test all(typeof.(rv0) == typeof.(sima))
    @test sima !== rv0
    # The expected behaviour (in some algorithms) is to drop the units here:
    simb = @inferred(similar(r0, Int64))
    @test all(typeof.(simb) .== Int64)
end

@testset "Inferrable similar vector and ArrayPartition mixed units" begin
    r0 = [1.0km, -2km, 3m/s, 4m/s]
    simi = @inferred(similar(r0))
    @test all(typeof.(r0) == typeof.(simi))
    @test simi !== r0
    rv0 = @inferred ArrayPartition(r0)
    sima = @inferred similar(rv0)
    @test all(typeof.(rv0) == typeof.(sima))
    @test sima !== rv0
    simb = @inferred(similar(r0, Int64))
    @test all(typeof.(simb) .== Int64)
end

@testset "Inferrable UNITLESS_ABS2 Unitless ArrayPartition" begin
    r0ul = [1131.340, -2282.343, 6672.423]
    v0ul = [-5.64305, 4.30333, 2.42879]
    rv0ul = ArrayPartition(r0ul,v0ul)
    @test @inferred(UNITLESS_ABS2(rv0ul)) == 5.101030471786125e7
end

@testset "Inferrable UNITLESS_ABS2  ArrayPartition mixed units" begin
    r0 = [1.0km, 2.0km, 3.0m/s, 4.0m/s]
    v0 = [1.0km/s, 2.0km/s, 3.0m/s², 4m/s²]
    rv0 = ArrayPartition(r0, v0)
    @test @inferred(UNITLESS_ABS2(1.0km)) == 1.0
    @test @inferred(UNITLESS_ABS2(r0)) == [1.0, 4.0, 9.0, 16.0]
    @test @inferred(UNITLESS_ABS2(rv0)) == 2 .* [1.0, 4.0, 9.0, 16.0]
end


text:
using Documenter, Schrodinger, SparseArrays, PyPlot

DocMeta.setdocmeta!(Schrodinger, :DocTestSetup, :(using Schrodinger, LinearAlgebra, SparseArrays); recursive=true)

makedocs(
    sitename = "Schrodinger.jl",
    authors = "Jérémy Béjanin.",
    modules = [Schrodinger],
    format = Documenter.HTML(prettyurls = false),
    linkcheck = true,
    clean = false,
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "man/gettingstarted.md",
            "man/quantumobjects.md",
            "man/working.md",
            "man/dynamics.md",
        ],
        "Examples" => [
            "examples/driven_TLS.md",
            "examples/DRAG.md",
        ],
        "API" => [
            "api/quobj.md",
            "api/states.md",
            "api/operators.md",
            "api/functions.md",
        ],
    ]
)

deploydocs(
    repo = "github.com/jebej/Schrodinger.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
module MultiChannelFactorAnalysis

import Base.+

export MultiChannelLayout, MultiChannelData, MultiChannelFactors, MultiChannelFactorLayout,  MCFMParams, MCFMHistory, MCFMFit
export fit, extract_factors, predict
export stack, unstack, extract_diagonal_blocks, extract_horizontal_blocks
export periodic_factors, classical_factors,indep_gauss_errors, simulate
export data_fromlong, to_long

using Statistics
using Distributions
using LinearAlgebra
using DataFrames
using Interpolations

include("structs.jl")
include("utils.jl")
include("MCFM.jl")
include("MCFMSim.jl")

end # module


text:
using BinaryBuilder

name = "SCALE_MAMBA"
version = v"1.5"

# Collection of sources required to build SuiteSparse
sources = [
    "https://github.com/KULeuven-COSIC/SCALE-MAMBA.git" =>
    "d7c960afd0a9776f04e15a5653caf300dd42f20a",
    "./bundled",
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/SCALE-MAMBA/

# Apply patch to fix multiple definitions error
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/fix_multiple_definitions.patch

cat > CONFIG.mine <<EOF
ROOT = $(pwd)
OSSL = ${prefix}

# Take default values from typical CONFIG
MAX_MOD = 7
FLAGS = -DMAX_MOD_SZ=\$(MAX_MOD)
OPT = -O3
EOF

if [[ ${target} == *linux* ]]; then
    echo "LDFLAGS = -lrt" >> CONFIG.mine
fi

# Build executables
make -j${nproc} CC="${CXX}" -C src

# Install executables into ${prefix}
mkdir -p ${prefix}/bin
cp Player.x${exe} ${prefix}/bin/
cp Setup.x${exe} ${prefix}/bin/
"""

# Only x86_64, no FreeBSD or windows, and no musl
platforms = [p for p in supported_platforms() if arch(p) == "x86_64" && !Sys.isfreebsd(p) && !Sys.iswindows(p) && libc(p) != "musl"]

# Build with GCC 6 at least, to dodge C++ problems
platforms = BinaryBuilder.replace_gcc_version.(platforms, :gcc6)

# The products that we will ensure are always built
products = [
    ExecutableProduct("Player.x", :playerx),
    ExecutableProduct("Setup.x", :setupx),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "OpenSSL_jll",
    "MPIR_jll",
    "cryptopp_jll",
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using MPIReco

@testset "motion compensation reconstruction" begin

  ################
  ## Parameters ##
  ################
  alpha = 1.2

  # Choose the peak number according to the expected approximate value. Be careful, if there are peaks not corresponding to object motion, division by choosePeak can be false.
  # For experiments with different periodic motions it is used to pick one motion frequency (then don't divide by choosePeak)

  choosePeak = 4

  windowType = 1  # 1: Hann, 2: FT1A05, 3: Rectangle

  # Measurement data
  datadirMeas = "./data/mdf/"
  bMeas = MPIFile(datadirMeas*"measFast.mdf") # high frequency
  bBG = MPIFile(datadirMeas*"measBG.mdf") # background measurement

  # System matrices
  datadirSF = "./data/mdf/"
  SFall = ["SF1Small.mdf","SF2Small.mdf","SF3Small.mdf","SF4Small.mdf"]
  bSF = MultiMPIFile(datadirSF.*SFall)

  # Background frames
  bgFrames = [1:200, 201:400, 401:600 ,601:800]

  # Selected frames
  recoFrame = 1

  ####################
  ## Reconstruction ##
  ####################

  # Reconstruction parameters
  lambda = 0.01
  iterations = 1

  c = reconstruction(bSF, bMeas, periodicMotionCorrection = true,
                                bEmpty = bBG, bgFrames = bgFrames, #background measurement
                                alpha = alpha, choosePeak = choosePeak, recoFrame = recoFrame,
                                windowType=windowType, higherHarmonic=choosePeak,
                                lambda=lambda, iterations=iterations, # reconstruction parameter
                                SNRThresh=10, minFreq=80e3, recChannels=[1,2,3]) #frequency selection

  exportImage("./img/MotionComp.png", maximum(arraydata(c[1,:,:,:,1]),dims=3)[:,:,1])

end


text:
module MappedArrays

using Base: @propagate_inbounds

export mappedarray

abstract AbstractMappedArray{T,N} <: AbstractArray{T,N}

immutable ReadonlyMappedArray{T,N,A<:AbstractArray,F} <: AbstractMappedArray{T,N}
    f::F
    data::A
end
immutable MappedArray{T,N,A<:AbstractArray,F,Finv} <: AbstractMappedArray{T,N}
    f::F
    finv::Finv
    data::A
end

"""
     mappedarray(f, A)

creates a view of the array `A` that applies `f` to every element of
`A`. The view is read-only (you can get values but not set them).
"""
mappedarray{T,N}(f, data::AbstractArray{T,N}) = ReadonlyMappedArray{typeof(f(one(T))),N,typeof(data),typeof(f)}(f, data)

"""
    mappedarray((f, finv), A)

creates a view of the array `A` that applies `f` to every element of
`A`. The inverse function, `finv`, allows one to also set values of
the view and, correspondingly, the values in `A`.
"""
function mappedarray{T,N}(f_finv::Tuple{Any,Any}, data::AbstractArray{T,N})
    f, finv = f_finv
    MappedArray{typeof(f(one(T))),N,typeof(data),typeof(f),typeof(finv)}(f, finv, data)
end

Base.parent(A::AbstractMappedArray) = A.data
Base.size(A::AbstractMappedArray) = size(A.data)
parenttype{T,N,A,F}(::Type{ReadonlyMappedArray{T,N,A,F}}) = A
parenttype{T,N,A,F,Finv}(::Type{MappedArray{T,N,A,F,Finv}}) = A
Base.linearindexing{MA<:AbstractMappedArray}(::Type{MA}) = Base.linearindexing(parenttype(MA))

@propagate_inbounds Base.getindex(A::AbstractMappedArray, i::Int...) = A.f(A.data[i...])
@propagate_inbounds Base.setindex!(A::MappedArray, val, i::Int...) = A.data[i...] = A.finv(val)

end # module


text:
@testset "supplier.jl" begin
    sup = Supplier(100, 10)
    @test state(sup) == []
    @test state_size(sup) == 0
    dest = []
    InventoryModels.pull!(sup, 13, dest)
    @test state(sup) == []
    @test only(sup.pull_orders) == (dest => 13)
    InventoryModels.dispatch!(sup)
    @test InventoryModels.reward!(sup) == -13*10-100
    @test dest == [13, sup.leadtime]

    sup2 = Supplier(FixedLinearOrderCost(100,10), capacity = 5)
    InventoryModels.pull!(sup2, 13, dest)
    @test only(sup2.pull_orders) == (dest => 5)
    @test state(sup2) == []
    InventoryModels.dispatch!(sup2)
    @test InventoryModels.reward!(sup2) == -5*10-100
    @test dest == [13, sup.leadtime, 5, sup2.leadtime]
    @test state(sup2) == []
    reset!(sup)

    sup = Supplier(100, 10, leadtime = LeadTime(1, 3, 10))
    @test state(sup) == [10,10,10]
    @test state_size(sup) == 3
    dest = []
    InventoryModels.pull!(sup, 13, dest)
    @test only(sup.pull_orders) == (dest => 13)
    InventoryModels.dispatch!(sup)
    @test state(sup) == [10,10,13]
    @test only(sup.leadtime.pull_orders) == (dest => 13)
    @test InventoryModels.reward!(sup) == -13*10-100 - 33
    @test dest == [10, sup.leadtime]
    @test InventoryModels.inventory_position(sup) == 33
end


text:
function solve_poisson_cg( LF::LF3dGrid, rho::Array{Float64,1}, NiterMax::Int64;
                           verbose=false, TOL=5.e-10 )
    Npoints = LF.Nx * LF.Ny * LF.Nz
    #
    phi = zeros( Float64, Npoints ) # XXX or use some starting guess
    #
    r = zeros( Float64, Npoints )
    p = zeros( Float64, Npoints )
    #
    nabla2_phi = apply_Laplacian( LF, phi )
    for ip = 1:Npoints
        r[ip] = rho[ip] - nabla2_phi[ip]
        p[ip] = r[ip]
    end

    rsold = dot( r, r )
    for iter = 1 : NiterMax
        #
        nabla2_phi = apply_Laplacian( LF, p )
        #
        alpha = rsold/dot( p, nabla2_phi )
        #
        phi = phi + alpha * p
        r = r - alpha * nabla2_phi
        #
        rsnew = dot( r, r )
        deltars = rsold - rsnew # used ?
        if verbose
          @printf("%8d %20.10f\n", iter, sqrt(rsnew))
        end
        #
        if sqrt(rsnew) < TOL
          if verbose
            @printf("#Convergence achieved in solve_poison_cg: N, iter: %d %d\n", Npoints, iter)
          end
          break
        end
        #
        p = r + (rsnew/rsold) * p
        #
        rsold = rsnew
    end
    #
    return phi
    #
end # of function

function solve_poisson_cg( Lmat::SparseMatrixCSC{Float64,Int64},
                           rho::Array{Float64,1}, NiterMax::Int64;
                           verbose=false, TOL=5.e-10 )
    #
    Npoints = size(rho)[1]
    phi = zeros( Float64, Npoints ) # XXX or use some starting guess
    #
    r = zeros( Float64, Npoints )
    p = zeros( Float64, Npoints )
    #
    nabla2_phi = Lmat*phi
    for ip = 1 : Npoints
        r[ip] = rho[ip] - nabla2_phi[ip]
        p[ip] = r[ip]
    end
    rsold = dot( r, r )

    for iter = 1 : NiterMax
        #
        nabla2_phi = Lmat*p
        #
        alpha = rsold/dot( p, nabla2_phi )
        #
        phi = phi + alpha * p
        r = r - alpha * nabla2_phi
        #
        rsnew = dot( r, r )
        deltars = rsold - rsnew
        if verbose
            @printf("%8d %20.10f\n", iter, sqrt(rsnew))
        end
        #
        if sqrt(rsnew) < TOL
            if verbose
                @printf("Convergence achieved in solve_poison_cg: N, iter: %d %d\n", Npoints, iter)
            end
            break
        end
        #
        p = r + (rsnew/rsold) * p
        #
        rsold = rsnew
    end
    #
    return phi
    #
end # of function