text:
using CompressedStacks

#import CompressedStacks.stack_test
import CompressedStacks.CompressedStack
#import CompressedStacks.access
#import CompressedStacks.run!

#stack_test(70,3)

#original benchmark function
#function benchmark()
    # Any setup code goes here.

    # Run once, to force compilation.
#    srand(666)
#    @time code_to_profile()

    # Run a second time, with profiling.
#    println("\n\n======================= Second run:")
#    srand(666)
#    Profile.init(delay=0.01)
#    clear_malloc_data()
#    @profile @time code_to_profile()

    # Write profile results to profile.bin.
#    r = Profile.retrieve()
#    f = open("profile.bin", "w")
#    serialize(f, r)
#    close(f)
#end

function printsum(a)
    # summary generates a summary of an object
    println(summary(a), ": ", repr(a))
end

function condition_push(stack::CompressedStack, elt::Int)

return true

end

function condition_pop(stack::CompressedStack, elt::Int)

  #println("La lio? ",isempty(stack))
  # This sort of reads the next one?
    if !isempty(stack)
      popRatio= parse(ARGS[3])
      min=parse(ARGS[4])
      max=parse(ARGS[5])
      range=max-min
      topElement=top(stack)

    #flush(STDOUT)
    #  println("should I POP ",aux[1]," ? : ",(parse(aux[1]) < (min+popRatio*range)))
    #flush(STDOUT)
        if (elt==(max+1))
          #println("POPPOP")
          return true
        else
          return ( (elt+topElement) < (min+2*popRatio*range))
        end
    else
        print(stack)
        println("IN this case the compressed stack was empty so I should NOT pop, unless I am comparing the run! of two stacks")
        return false
    end

end

function action_push(stack::CompressedStack{Int,Int}, elt::Int)
  stack.context = Nullable(stack.index)
  # print(stack)
end

function action_pop(stack::CompressedStack, elt::Int)
#  println("Pop element : $elt")
end

#now do the same thing again for the normal stack

function condition_push2(stack::NormalStack, elt::Int)
#  aux = elt

#  popRatio= parse(ARGS[3])
#  min=parse(ARGS[4])
#  max=parse(ARGS[5])
#  range=max-min
#  return (aux > (min+popRatio*range))
  return true
end

function condition_pop2(stack::NormalStack, elt::Int)
#println("La lio? ",isempty(stack))
# This sort of reads the next one?
  if !isempty(stack)
    popRatio= parse(ARGS[3])
    min=parse(ARGS[4])
    max=parse(ARGS[5])
    range=max-min
    topElement=top(stack)

  #flush(STDOUT)
  #  println("should I POP ",aux[1]," ? : ",(parse(aux[1]) < (min+popRatio*range)))
  #flush(STDOUT)
      if (elt==(max+1))
          #println("POPPOP")
          return true
      else
        return ( (elt+topElement) < (min+2*popRatio*range))
      end

  else
      print(stack)
      println("IN this case the normal stack was empty so I should NOT pop, unless I am comparing the run! of two stacks")
      return false
  end

end

function action_push2(stack::NormalStack{Int,Int}, elt::Int)
  stack.context = Nullable(stack.index)
  # print(stack)
end

function action_pop2(stack::NormalStack, elt::Int)
#  println("Pop element : $elt")
end

function randomPushPopTestCompressed(name,size,space)

  println("randomPushPopTest with compressed stack, with Parameters size:",size," space: ",space," name ",name)

  context_type = Int
  data_type = Int

  #ns = NormalStack(name, action_pop2, action_push2, condition_pop2,
  #condition_push2, context_type, data_type)

  #cs = CompressedStack(name, action_pop, action_push, condition_pop,
  #condition_push, context_type, data_type)

  stack = CompressedStack(name, action_pop, action_push, condition_pop, condition_push, context_type, data_type,size,space)
 normalstack = NormalStack(name, action_pop2, action_push2, condition_pop2, condition_push2, context_type, data_type)

#println("Running ")
#try run!(stack,normalstack)
#try  run!(stack)
#try run!(normalstack)
#catch Exception e
#  println("Yagos code caught this exception ",e)
#  println("The stack was ")
#  print(stack)
#end

# seekend(stack.input)
# seekend(normalstack.input)
# println(position(stack.input)==position(normalstack.input))

run!(stack,normalstack)
#run!(stack)

  stack = 0

end

#function randomPushPopTest(repetitions,popRatio,min,max)
function randomPushPopTest(name)

# First let us make a normal stack
#q=Int[]
#println("randomPushPopTest, Parameters ",repetitions," ",popRatio," ",min," ",max," ")
#range=max-min
# Careful, as the rand is of integers this is not really working as a random choice with popRatio parameter
#for i = 1:repetitions
#  (x=rand(min:max)) > (min+popRatio*range) ? push!(q,x) : isempty(q) ? true : pop!(q)
  #println(i)
#end

println("randomPushPopTest with normal stack, with Parameters name:",name)

context_type = Int
data_type = Int

#ns = NormalStack(name, action_pop2, action_push2, condition_pop2,
#condition_push2, context_type, data_type)

#cs = CompressedStack(name, action_pop, action_push, condition_pop,
#condition_push, context_type, data_type)

#stack = CompressedStack(name, action_pop, action_push, condition_pop, condition_push, context_type, data_type,size,space)
normalstack = NormalStack(name, action_pop2, action_push2, condition_pop2, condition_push2, context_type, data_type)

run!(normalstack)
#try run!(normalstack)
#catch Exception e
#  println("normal stack caught this exception ",e)
#  println("The stack was ")
#  print(normalstack)
#end

normalstack = 0

end

#printsum(q)
if (isempty(ARGS) || length(ARGS)<8) println("WRONG PARAMETERS! NEED 5 (repetitions popRatio min max n p testType(0 normal 1 compressed) )")
else

  IOFileName=ARGS[1]
  size=parse(ARGS[2])
  min=parse(ARGS[4])
  max=parse(ARGS[5])
  n=parse(ARGS[6])
  p=parse(ARGS[7])

   println("Testing : FileName ",IOFileName,"\n repetitions: ",size," popRatio  ",ARGS[3]," min ",min,"  max ",max,"  n ",n," p ",p," testType(0 normal 1 compressed) ",ARGS[8])

    # Run once, to force compilation.
  #println("======================= First run:")
  #srand(666)
  #@time parse(ARGS[8])==0 ? randomPushPopTest(10,0.1,1,5) : randomPushPopTestCompressed("/home/yago/.julia/v0.4/CompressedStacks/ioexample/input1")

  # Run a second time, with profiling.
  #println("\n\n======================= Second run:")
  srand(666)

  Profile.init(delay=0.01)
  Profile.clear()

  Profile.clear_malloc_data()
  gc()
  #@profile @time parse(ARGS[8])==0 ? randomPushPopTest(parse(ARGS[2]),parse(ARGS[3]),parse(ARGS[4]),parse(ARGS[5])) : randomPushPopTestCompressed(IOFileName)
#parse(ARGS[8])==0 ? randomPushPopTest(parse(ARGS[2]),parse(ARGS[3]),parse(ARGS[4]),parse(ARGS[5])) : randomPushPopTestCompressed(IOFileName,n,p)
parse(ARGS[8])==0 ? randomPushPopTest( IOFileName ) : randomPushPopTestCompressed(IOFileName,n,p)

  # Write profile results to profile.bin.
#  r = Profile.retrieve()
#  f = open("profile.bin", "w")
#  serialize(f, r)
#  close(f)

end


text:
# run this file using julia --project -t4 ch14_server.jl command

using Genie
using Statistics
using ThreadsX

function payoff_asian_sample(T, X0, K, r, s, m)::Float64
    X = X0
    sumX = X
    d = T / m
    for i in 1:m
        X *= exp((r - s^2 / 2) * d + s * sqrt(d) * randn())
        sumX += X
    end
    Y = sumX / (m + 1)
    return exp(-r * T) * max(Y - K, 0)
end

function asian_value(T, X0, K, r, s, m, max_time)
    result = Float64[]
    start_time = time()
    while time() - start_time < max_time
        append!(result,
                ThreadsX.map(_ -> payoff_asian_sample(T, X0, K, r, s, m),
                             1:10_000))
    end
    n = length(result)
    mv = mean(result)
    sdv = std(result)
    lo95 = mv - 1.96 * sdv / sqrt(n)
    hi95 = mv + 1.96 * sdv / sqrt(n)
    zero = mean(==(0), result)
    return (; n, mv, lo95, hi95, zero)
end

Genie.config.run_as_server = true

Genie.Router.route("/", method=POST) do
  message = Genie.Requests.jsonpayload()
  return try
      K = float(message["K"])
      max_time = float(message["max_time"])
      value = asian_value(1.0, 50.0, K, 0.05, 0.3, 200, max_time)
      Genie.Renderer.Json.json((status="OK", value=value))
  catch
      Genie.Renderer.Json.json((status="ERROR", value=""))
  end
end

Genie.startup()

ext: jl

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
using InferenceUtilities
using Base.Test

# tests copied from pull request
# https://github.com/JuliaLang/julia/pull/23426

# test @inferred and @isinferred
function uninferrable_function(i)
    q = [1, "1"]
    return q[i]
end

@test_throws ErrorException @inferred(uninferrable_function(1))
@test @inferred(identity(1)) == 1
@test !@isinferred(uninferrable_function(1))
@test @isinferred(identity(1))

# Ensure @inferred and @isinferred only evaluate the arguments once
inferred_test_global = 0
function inferred_test_function()
    global inferred_test_global
    inferred_test_global += 1
    true
end
@test @inferred inferred_test_function()
@test inferred_test_global == 1

inferred_test_global = 0
@test @isinferred inferred_test_function()
@test inferred_test_global == 1

# Test that @inferred and @isinferred work with A[i] expressions
@test @inferred((1:3)[2]) == 2
@test @isinferred((1:3)[2])
struct SillyArray <: AbstractArray{Float64,1} end
Base.getindex(a::SillyArray, i) = rand() > 0.5 ? 0 : false
test_result = @test_throws ErrorException @inferred(SillyArray()[2])
@test contains(test_result.value.msg, "Bool")
@test !@isinferred(SillyArray()[2])

# Issue #14928
# Make sure abstract error type works.
@test_throws Exception error("")

# Issue #17105
# @inferred and @isinferred with kwargs
function inferrable_kwtest(x; y=1)
    2x
end
function uninferrable_kwtest(x; y=1)
    2x+y
end
@test @inferred(inferrable_kwtest(1)) == 2
@test @inferred(inferrable_kwtest(1; y=1)) == 2
@test @inferred(uninferrable_kwtest(1)) == 3
@test_throws ErrorException @inferred(uninferrable_kwtest(1; y=2)) == 2
@test @isinferred(inferrable_kwtest(1))
@test @isinferred(inferrable_kwtest(1; y=1))
@test @isinferred(uninferrable_kwtest(1))
@test !@isinferred(uninferrable_kwtest(1; y=2))


text:
module GaussianMixturesController
using StatsBase, LinearAlgebra
using ..Stats, ..GaussianMixture
using JSON

export call

"""
call(message)

Computes gaussian mixtures from records.

# Arguments
- `message::Dict`: the arguments dict
  - `"records"`: the array of array of rows
  - `"components"`: component number (limitation: less than or equal to the row size)

# Return
the json string of result object

# Examples
```jldoctest
julia> GaussianMixturesController.call(Dict("components" => 2, "records" => [[1.0, 2.0, 0.0], [2.0, 1.0, 2.0], [3.0, 0.0, 9.0]]))
"{\"means\":[[2.499966642114253,0.9999999999999978],[0.5000333578858548,1.9999999999999956],[0.9999999999999991,8.99999999999998]],\"weights\":[[0.6666814927233754,0.3333185072766262],[0.6666814927234248,0.33331850727657647],[0.6666666666666674,0.3333333333333341]],\"stds\":[[0.5000454740003658,0.001000000001110223],[0.5000454740005074,0.0010000000044408921],[1.0000004999998753,0.0010000000852651246]]}"
```
"""
function call(message)
  records = hcat(convert.(Array{Any,1}, message["records"])...)
  components = message["components"]
  means, stds, weights, mins, maxs = [], [], [], [], []
  for i in 1:size(records, 1)
    values = convert(Array{Float64, 1}, view(records, i, :))
    best_gmm = GaussianMixture.fit(values, components)

    push!(means, best_gmm.means_[:, 2])
    covars = best_gmm.covariances_[:, 2]
    push!(stds, sqrt.(covars))
    push!(weights, best_gmm.weights_)
  end

  Dict(
       :means => means,
       :stds => stds,
       :weights => weights,
      ) |> JSON.json
end
end


text:
# based on http://int8.io/automatic-differentiation-machine-learning-julia/

using MNIST # if not installed try Pkg.add("MNIST")
using AutoDiffSource # if not installed try Pkg.add("AutoDiffSource")
using PyPlot # if not installed try Pkg.add("PyPlot")

sigmoid(x) = 1 / (1 + exp(x))

function autoencoder(We1, We2, Wd, b1, b2, input)
    firstLayer = sigmoid.(We1 * input + b1)
    encodedInput = sigmoid.(We2 * firstLayer + b2)
    reconstructedInput = sigmoid.(Wd * encodedInput)
    return reconstructedInput
end

@δ function autoencoderError(We1, We2, Wd, b1, b2, input)
    reconstructedInput = autoencoder(We1, We2, Wd, b1, b2, input)
    return sum((input - reconstructedInput).^2)
end

@assert checkdiff(autoencoderError, δautoencoderError, randn(3,3), randn(3,3), rand(3,3), randn(3), randn(3), randn(3))

function initializeNetworkParams(inputSize, layer1Size, layer2Size)
    We1 =  1f-1 * randn(Float32, layer1Size, inputSize)
    b1 = zeros(Float32, layer1Size)
    We2 =  1f-1 * randn(Float32, layer2Size, layer1Size)
    b2 = zeros(Float32, layer2Size)
    Wd = 1f-1 * randn(Float32, inputSize, layer2Size)
    return (We1, We2, b1, b2, Wd)
end

function show_digits(testing, We1, We2, b1, b2, Wd)
    clf()
    total_error = 0
    for l = 1:12
        input = testing[:, rand(1:size(testing, 2))]
        reconstructedInput = autoencoder(We1, We2, Wd, b1, b2, input)
        subplot(4, 6, l*2-1)
        pcolor(rotl90(reshape(input, 28, 28)'); cmap="Greys")
        subplot(4, 6, l*2)
        pcolor(rotl90(reshape(reconstructedInput, 28, 28)'); cmap="Greys")
        total_error += sum((input - reconstructedInput).^2)
    end
    tight_layout()
    total_error / 12
end

function trainAutoencoder(epochs, training, testing, We1, We2, b1, b2, Wd, alpha)
    for k in 1:epochs
        total_error = 0.
        for i in 1:size(training, 2)
            input = training[:,i]
            val, ∇autoencoderError = δautoencoderError_const6(We1, We2, Wd, b1, b2, input)
            total_error += val
            if mod(i, 10_000) == 0
                test_error = show_digits(testing, We1, We2, b1, b2, Wd)
                @printf("epoch=%d iter=%d train_error=%.2f test_error=%.2f\n", k, i, total_error/10_000, test_error)
                total_error = 0.
            end
            ∂We1, ∂We2, ∂Wd, ∂b1, ∂b2 = ∇autoencoderError()
            We1 -= alpha * ∂We1
            We2 -= alpha * ∂We2
            Wd  -= alpha * ∂Wd
            b1  -= alpha * ∂b1
            b2  -= alpha * ∂b2
        end
    end
    return (We1, We2, b1, b2, Wd)
end

# read input MNIST data
training = convert(Matrix{Float32}, MNIST.traindata()[1] / 255)
testing = convert(Matrix{Float32}, MNIST.testdata()[1] / 255)

# 784 -> 300 -> 100 -> 784 with weights normally distributed (with small variance)
We1, We2, b1, b2, Wd = initializeNetworkParams(784, 300, 100)

# 4 epochs with alpha = 0.02
@time We1, We2, b1, b2, Wd = trainAutoencoder(4, training, testing, We1, We2, b1, b2, Wd, 2f-2)


text:
export Predicate
export all
export exists

"""
    Represents a predicate.

A proposition-valued function of one or more variables.

fields: value and x.
"""
struct Predicate <: Logic
    value::Logic
    x::Tuple
    Predicate(p::Logic, x...) = new(p, tuple(x...))
    Predicate(P::Predicate, R::Predicate) = new(P.value * R.value, (P.x..., R.x...))
end

"""
    all(x, P)

Construct a proposition with the given Predicate `P` and a single variable `x`.
"""
all(x, P::Predicate) = value(P)

"""
    exists(x, P)

Construct a proposition with the given Predicate `P` and a single variable `x`.
"""
exists(x, P::Predicate) = !all(x, !P)


text:
## run app.jl then define the following function
myhook = true
function button_test(store)
    while myhook
        if get_state(store).buttons["basic_button"].button.is_triggered
            @warn "This triggers $(@__FILE__):$(@__LINE__)."
        end
        yield()
    end
end

## then call the function asyncly
t = @async button_test(store)

# click the button to see what happens

## then set `myhook` to false and the function will exit
myhook = false
@show t

## update the function and turn on `myhook` again
myhook = true
t = @async button_test(store)

## this is basically how to interactively connect/disconnect functions from the UI
# there is no callback functions or signals involved, all we need is just a well-managed
# global state along with naive Julia tasks.


text:
"""
    std_error(x[; method])

Estimate the standard error of the mean of the given time series.

The keyword `method` can be used to choose one out of the following methods:

* Logarithmic binning: `:log` (default)
* Full binning: `:full`
* Jackknife resampling: `:jackknife`
"""
function std_error(x::AbstractVector{T}; method::Symbol=:log) where T
    if method == :log
        B = LogBinner(zero(x[1]), capacity=length(x))
        append!(B, x)
        return std_error(B)
    elseif method == :full
        return std_error(FullBinner(x))
    elseif method == :jackknife
        T <: Number || error("Jackknife doesn't support non-number type time series.")
        return Jackknife.std_error(identity, x)
    elseif method == :error_propagator
        B = ErrorPropagator(zero(x[1]), capacity=length(x))
        append!(B, x)
        return std_error(B, 1)
    else
        throw(ArgumentError(
            "Keyword `method` must be either `:log`, `:full`, `:jackknife` or " *
            "`:error_propagator`. Got `$(method)`."
        ))
    end
end

"""
    tau(x[; method])

Estimate the autocorrelation time of the given time series.

The keyword `method` can be used to choose one out of the following methods:

* Logarithmic binning: `:log` (default)
* Full binning: `:full`
"""
function tau(x::AbstractVector{T}; method::Symbol=:log) where T
    if method == :log
        B = LogBinner(zero(x[1]), capacity=length(x))
        append!(B, x)
        return tau(B)
    elseif method == :full
        return tau(FullBinner(x))
    else
        throw(ArgumentError("Keyword `method` must be either `:log` or `:full`. Got `$(method)`."))
    end
end


text:
push!(LOAD_PATH, pwd())

module TestModule

using A
using B
using Test
using JLD2

x =  AType(1)

S = BType(x)

@testset "save and load" begin
	T = deepcopy(S)
	@test T == S
	@save "test_out.jld" T
	T = nothing
	@load "test_out.jld" T
	@test T == S
end

end


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).
using Test

using SolidStateDetectors
using SolidStateDetectors: SSDFloat

using Unitful

T = Float32

@testset "Comparison to analytic solutions" begin
    include("comparison_to_analytic_solutions.jl")
end

@testset "Test real detectors" begin
    @testset "Simulate example detector: Inverted Coax" begin
        sim = Simulation(SSD_examples[:InvertedCoax])
        simulate!(sim, max_refinements = 1, verbose = true)
        evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 40e-3 )]))
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.waveforms)
            signalsum += abs(evt.waveforms[i].value[end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end
    @testset "Simulate example detector: Coax" begin
        sim = Simulation(SSD_examples[:Coax])
        SolidStateDetectors.apply_initial_state!(sim, ElectricPotential)
        # simulate!(sim, max_refinements = 0, verbose = true)
        # evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(30), 12e-3 )]))
        # simulate!(evt, sim)
        # signalsum = T(0)
        # for i in 1:length(evt.waveforms)
        #     signalsum += abs(evt.waveforms[i].value[end])
        # end
        # @test isapprox( signalsum, T(2), atol = 5e-2 )
    end
    @testset "Simulate example detector: BEGe" begin
        sim = Simulation(SSD_examples[:BEGe])
        simulate!(sim, max_refinements = 0, verbose = true)
        evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 20e-3 )]))
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.waveforms)
            signalsum += abs(evt.waveforms[i].value[end])
        end
        @test isapprox( signalsum, T(2), atol = 1e-2 )
    end
    @testset "Simulate example detector: CGD" begin
        sim = Simulation(SSD_examples[:CGD])
        SolidStateDetectors.apply_initial_state!(sim, ElectricPotential)
        # simulate!(sim, max_refinements = 0, verbose = true)
        # evt = Event([CartesianPoint{T}(5e-3, 5e-3, 5e-3)])
        # simulate!(evt, sim)
        # signalsum = T(0)
        # for i in 1:length(evt.waveforms)
        #     signalsum += abs(evt.waveforms[i].value[end])
        # end
        # @test isapprox( signalsum, T(2), atol = 1e-2 )
    end
    @testset "Simulate example detector: Spherical" begin
        sim = Simulation(SSD_examples[:Spherical])
        SolidStateDetectors.apply_initial_state!(sim, ElectricPotential)
        # simulate!(sim, max_refinements = 1, verbose = true)
        # evt = Event([CartesianPoint{T}(0,0,0)])
        # simulate!(evt, sim)
        # signalsum = T(0)
        # for i in 1:length(evt.waveforms)
        #     signalsum += abs(evt.waveforms[i].value[end])
        # end
        # @test isapprox( signalsum, T(2), atol = 1e-2 )
    end
    @testset "Simulate example detector: SigGen PPC" begin
        sim = Simulation(SSD_examples[:SigGen])
        simulate!(sim, max_refinements = 0, verbose = true)
        evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 40e-3 )]))
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.waveforms)
            signalsum += abs(evt.waveforms[i].value[end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end
end


text:
# ZJ: currenlty there is no support for distributed array and given the speed below it's not ready
using StatsBase

@time addprocs()
@everywhere using FastGroupBy

@time did6 = SharedArray(id6)
@time dv1 = SharedArray(v1)

function sumby_sa(by, val)
  T = eltype(by)
  S = eltype(val)

  l = length(by)
  a = sort(collect(Set([1:Int(round((l/nworkers()))):l...,l])))

  if length(a) > nworkers() + 1
    a_mid = a[2:(end-1)]
    a = vcat(1, sort(setdiff(a_mid, sample(a_mid, length(a) - nwokers() - 1, repalce = false))), l)
  end
  b = a[2:end]

  res = [@spawnat k sumby(by[a[k]:b[k]],val[a[k]:b[k]]) for k = 1:nworkers()]

  @time res1 = fetch.(res)
  return res1

  # algorithms to collate all dicts
  fnl_res = fetch(res[1])
  szero = zero(S)
  @time for i = 2:length(res)
    next_res = fetch(res[i])
    for k = keys(next_res)
      fnl_res[k] = get(fnl_res, k, szero) + next_res[k]
    end
  end
  fnl_res
end

@time sumby_sa(did6,dv1)


text:
using LoopVectorization, Test, ArrayInterface
using LoopVectorization: check_inputs!

# macros for generate loops whose body is not a block
macro gen_loop_issue395(ex)
  sym, ind = ex.args
  loop_body = :(ret[$ind] = $sym[$ind])
  loop = Expr(:for, :($ind = axes($sym, 1)), loop_body)
  return esc(:(@turbo $loop))
end
macro gen_single_loop(B, A)
  loop_body = :($B[i] = $A[i])
  loop = Expr(:for, :(i = indices(($B, $A), 1)), loop_body)
  return esc(:(@turbo $loop))
end
macro gen_nest_loop(C, A, B)
  loop_body = :($C[i, j] = $A[i] * $B[j])
  loop_head = Expr(:block, :(j = indices(($C, $B), (2, 1))), :(i = indices(($C, $A), 1)))
  loop = Expr(:for, loop_head, loop_body)
  return esc(:(@turbo $loop))
end
macro gen_A_mul_B(C, A, B)
  inner_body = :(Cji += $A[j, k] * $B[k, i])
  inner_loop = Expr(:for, :(k = indices(($A, $B), (2, 1))), inner_body)
  loop = :(
    for i in indices(($C, $B), 2), j in indices(($C, $A), 1)
      Cji = zero(eltype($C))
      $inner_loop
      $C[j, i] = Cji
    end
  )
  return esc(:(@turbo $loop))
end

@testset "check_block, #395" begin
  A = rand(4)
  B = rand(4)
  C = rand(4, 4)
  D = zeros(4)
  E = zeros(4, 4)
  F = zeros(4, 4)
  ret = zeros(4)
  @gen_single_loop D A
  @gen_loop_issue395 B[i]
  @gen_nest_loop E A B
  @gen_A_mul_B F C E
  @test D == A
  @test ret == B
  @test E == A * B'
  @test F == C * E
end

@testset "enumerate, #393" begin
  A = zeros(4)
  B = zeros(4)
  C = zeros(4, 4)
  D = zeros(4, 4)
  @turbo for (i, x) in enumerate(1:4)
    A[i] = x
  end
  @turbo for (i,) in enumerate(B)
    B[i] += 1
  end
  @turbo for (j, Aj) in enumerate(A), (i, Bi) in enumerate(B)
    C[i, j] = Aj * Bi
  end
  @turbo for (j, Bj) in enumerate(B)
    for (i, Ai) in enumerate(A)
      D[i, j] = Ai * Bj
    end
  end
  @test A == 1:4
  @test B == ones(4)
  @test A .* B' == C' == D
  @test_throws ArgumentError check_inputs!(:(
    for ix in enumerate(A)
      A[ix[1]] = ix[1] + ix[2]
    end
  ), Any[])
  @test_throws ArgumentError check_inputs!(:(
    for () in enumerate(A)
    end
  ), Any[])
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # The gradient regularized asymptotic correction.
# ## I. The ionization potenital issue
#
# Often with density functional theory one would like to do something interesting with the resulting wavefunction; however a number of issues exsist with conventional DFT computations.
#
# Let us first compute Ionization Potential at several different levels of the theory, and then compare the results to the NIST database (http://physics.nist.gov).

# +
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
using LinearAlgebra: dot
import Formatting: printfmt
include("ks_helper.jl")
import PyPlot: plt, matplotlib
matplotlib.use(backend="MacOSX")

build_superfunctional = nothing
if VersionNumber(psi4.__version__) >= v"1.3a1"
    build_superfunctional = psi4.driver.dft.build_superfunctional
else
    build_superfunctional = psi4.driver.dft_funcs.build_superfunctional
end

mol = psi4.geometry("""
Ne
symmetry c1
""")
options = (Dict("BASIS" => "aug-cc-pvdz",
                "NAT_ORBS" => true,
                "NUM_ROOTS" => 2,
                "QC_MODULE" => "DETCI",
                "DFT_SPHERICAL_POINTS" => 110,
                "DFT_RADIAL_POINTS" => 20))
psi4.set_options(options)

scf_e, scf_wfn = psi4.energy("SCF",  return_wfn=true)
scf_e, dft_wfn = psi4.energy("PBE0", return_wfn=true)
scf_e, ci_wfn =  psi4.energy("CISD", return_wfn=true)
# -

# Now that we have the wavefunctions we can compute the ionization potenital. For SCF theories this is the energy of the HOMO while for more advanced theories, it is the differente between the ground and excited states. We can obtain the IP for each level of theory in the following way:

# +
scf_ip = -1 * scf_wfn.epsilon_a().np[scf_wfn.nalpha()] * psi4.constants.hartree2ev
dft_ip = -1 * dft_wfn.epsilon_a().np[dft_wfn.nalpha()] * psi4.constants.hartree2ev
ci_ip = -1 * (psi4.core.variable("CI ROOT 0 TOTAL ENERGY") -
              psi4.core.variable("CI ROOT 1 TOTAL ENERGY")) * psi4.constants.hartree2ev

println("Neon Ionization Potential (eV)")
printfmt("SCF     {:10.4f}\n", scf_ip)
printfmt("PBE0    {:10.4f}\n", dft_ip)
printfmt("CISD    {:10.4f}\n", ci_ip)
printfmt("NIST    {:10.4f}\n", 21.5645)

# -

# ## I. Density shapes
# If we think carefully about the ionization potential and what this really means, it can be summed up susinctly as "How hard is it for the electron to leave the system" or "How deep is our potenital well"? A easy way to view the shape of the potenital for a given method is to inspect the density as a function of spatial extent. A weaker potential would mean that the density extends further than it should, and our IP would be too small. Whereas a stronger potenital would mean that the density is more concentrated and the IP is too large. As we only have a single atom we can project our density into a spatial grid and then sort the results as a function of distance from the nucleus.

# +
function one_electron_radial(matrices, Vpot)
    """
    For every matrix in matrices, computes the radial density using the incoming Vpot object.
    """

    results = [[] for x in matrices]
    Rlist = []

    points_func = Vpot.properties()[1]
    superfunc = Vpot.functional()

    xc_e = 0.0
    vv10_e = 0.0

    # First loop over the outer set of blocks
    for block in 1:Vpot.nblocks()

        # Obtain general grid information
        grid = Vpot.get_block(block-1)
        x = np.array(grid.x())
        y = np.array(grid.y())
        z = np.array(grid.z())
        npoints = size(z,1)

        # Get radial information
        R = @. (x^2 + y^2 + z^2)^0.5
        push!(Rlist, R)

        # Loop over and compute the density for each matrix
        for (num, mat) in enumerate(matrices)
            points_func.set_pointers(mat)
            points_func.compute_points(grid)
            push!(results[num], np.array(points_func.point_values()["RHO_A"])[1:npoints])
        end
    end

    # Stack R
    Rlist = vcat(Rlist...)

    # Get the R unique and sort
    #uidx = findfirst.(isequal.(unique(Rlist)), [Rlist])
    uidx = indexin(sort(unique(Rlist)), Rlist)
    Rlist = Rlist[uidx]

    # Julia does magic, it would be worthwhile to consider what this does exactly
    results = [vcat(x...)[uidx] for x in results]

    Rlist, results
end

# Compute the radial densities
sup = build_superfunctional("SVWN", true)[1]
Vpot = psi4.core.VBase.build(scf_wfn.basisset(), sup, "RV")
Vpot.initialize()

R, data = one_electron_radial([scf_wfn.Da(), dft_wfn.Da(), ci_wfn.Da()], Vpot);
# -

# Now that we have the radial densities, we build  a function that will plot the density as a ratio between
# the method and a benchmark. Here we will select the CISD densities as the benchmark and plot the SCF and DFT quantities against them.

function plot_scatter_ratio(dictionary, bench, x_values)
    fig1 = plt.figure()
    ax1 = fig1.add_subplot(111)
    for (k, v) in dictionary
        ax1.plot(x_values, v ./ bench, label=k)
    end

    ax1.legend(loc="upper left")
    ax1.set_xlim(0, 7)
    ax1.set_xlabel("R")
    ax1.set_ylabel(raw"ρ(r) / ρ(r)_bench")

    plt.plot([0, 7], [1.0, 1.0], "k-")
end

# Plotting this data we note that the SCF density falls off too rapidly, meaning that potential is too deep while the DFT density falls off far too slowly indicating the potenital is too shallow. For example many theories examine the overlap between densities. From this chart we can surmise issues that would arise with both DFT and HF theories.

data_dict = Dict("SCF" => data[1], "DFT" => data[2])
plot_scatter_ratio(data_dict, data[3], R)

# ## 2. The GRAC Correction
# A question then arises: what if we could correct the DFT potential asymptote so that it recovers the correct density and therefore the correction ionization potenital? Effectively, we want to accomplish:
#
# $$v_{\rm{xc}}({\bf r}) = -\frac{1}{r} + v_{\rm{xc}}(\infty).$$
#
# However, having a correct asympote would require knowledge of the exact $v_{\rm{xc}}$ potential. As we operate completely in the Density Funcitonal Approximation domain, this is not a possible task. Instead, the GRAC correction requires that we know our Ionization Potential and we can replace the complex DFT term with the following:
#
# $$v_{\rm{xc}}({\bf r}) = -\frac{1}{r} + (I_{\bf P} + \epsilon_N)$$
#
# Where $\epsilon_N$ is our HOMO energy. This is quite clever as we now only need a functional that correctly goes to the $\frac{1}{r}$ asymptotic. As no general functional accomplishes this task a special functional, "LB94", was built with exactly this property. The downside is that the characteristics of this functional are quite poor in the "bulk" region, or the region near the atom itself. To generalize this a switching function is introduced:
#
# $$
# \begin{align}
# v_{\rm{xc}}^{GRAC}(\rho({\bf r})) &= [1 - f[x({\bf r})]]v_{\rm{XC}}^{\rm bulk}+ f[x({\bf r})]v_{\rm{XC}}^{\rm Asym.}  \\ \\
# f[x({\bf r})] &= \frac{1}{1 + e^{-\alpha[x({\bf r}) - \beta]}}\\\\
# x({\bf r}) &= \frac{| \nabla \rho({\bf r}) | }{\rho^\frac{4}{3}({\bf r})}
# \end{align}
# $$
#
#
# We seamlessly switch between the bulk and asymptotic potenitals. In addition, such a formulation allows us to apply the GRAC correction to any underlying functional in a general way.
#

# +
# Build a LB 94 functional
lb94_func = psi4.core.SuperFunctional.blank()
lb94_func.set_max_points(500000)
lb94_func.set_deriv(1)

# Add member functionals
tmp_func = psi4.core.LibXCFunctional("XC_GGA_X_LB", true)
# We need to scale this for the amount of exact exchange is in the functional
tmp_func.set_alpha(0.75)
lb94_func.add_x_functional(tmp_func)

# Dont forget VWN correlation!
tmp_func = psi4.core.LibXCFunctional("XC_LDA_C_VWN", true)
lb94_func.add_c_functional(tmp_func)

# Allocate the functional
lb94_func.allocate()

grac_shift = (21.5645 - dft_ip) /  psi4.constants.hartree2ev

"""
Computes the XC energy and the V_xc Fock term
"""
function compute_V_GRAC(D, Vpot)

    grac_alpha = 0.5
    grac_beta = 40.0

    nbf = D.shape[1]
    Varr = zeros(nbf, nbf)

    points_func = Vpot.properties()[1]
    superfunc   = Vpot.functional()

    total_e = 0.0

    for x in 1:Vpot.nblocks()

        # Grid data
        grid = Vpot.get_block(x-1)
        w = np.array(grid.w())
        npoints = size(w,1)

        points_func.compute_points(grid)
        lpos = np.array(grid.functions_local_to_global()) .+ 1
        nfunctions = size(lpos,1)

        phi = np.array(points_func.basis_values()["PHI"])[1:npoints, 1:nfunctions]

        rho = np.array(points_func.point_values()["RHO_A"])[1:npoints]
        gamma = np.array(points_func.point_values()["GAMMA_AA"])[1:npoints]

        grac_x = @. (gamma^0.5) / (rho^(4.0/3.0))
        grac_fx = @. 1.0 / (1.0 + exp(-grac_alpha * (grac_x - grac_beta)))

        ret = superfunc.compute_functional(points_func.point_values(), -1)
        lbret = lb94_func.compute_functional(points_func.point_values(), -1)

        total_e += dot(np.array(ret["V"])[1:npoints], w)

        v_rho_a  = @. ( 1 - grac_fx ) * (np.array(ret["V_RHO_A"])[1:npoints] - grac_shift)
        v_rho_a += @. grac_fx * np.array(lbret["V_RHO_A"])[1:npoints]
        v_rho_a .*= 0.5 .* w

        # Vtmp[ab] = ϕ[pb] Vρ_a[p] ϕ[pa]
        Vtmp = phi' * (v_rho_a .* phi)

        # GGA
        rho_x = np.array(points_func.point_values()["RHO_AX"])[1:npoints]
        rho_y = np.array(points_func.point_values()["RHO_AY"])[1:npoints]
        rho_z = np.array(points_func.point_values()["RHO_AZ"])[1:npoints]

        phi_x = np.array(points_func.basis_values()["PHI_X"])[1:npoints, 1:nfunctions]
        phi_y = np.array(points_func.basis_values()["PHI_Y"])[1:npoints, 1:nfunctions]
        phi_z = np.array(points_func.basis_values()["PHI_Z"])[1:npoints, 1:nfunctions]

        v_gamma  = @. (1 - grac_fx) * np.array(ret["V_GAMMA_AA"])[1:npoints]
        v_gamma .*= 2.0 .* w
        # Vtmp[ab] += ϕ[pb] v_gamma[p] ∇ₓρ[p] ∇ₓϕ[pa] ∀ x,y,z
        Vtmp += phi' * (v_gamma .* rho_x .* phi_x)
        Vtmp += phi' * (v_gamma .* rho_y .* phi_y)
        Vtmp += phi' * (v_gamma .* rho_z .* phi_z)

        Varr[lpos, lpos] += Vtmp + Vtmp'

    end
    total_e, Varr
end

# Test our function the the handy ks solver.
grac_e, grac_data = ks_solver("PBE0", mol, options, compute_V_GRAC)
# -

# We can now check against the Psi4 reference implementation. Note that we set the GRAC shift on either side of the computation so that we do not contaminate other computations.

psi4.set_options(Dict("dft_grac_shift" => grac_shift))
println("Psi4 PBE0 energy ", psi4.energy("PBE0"))
psi4.set_options(Dict("dft_grac_shift" => 0.0))

# ## 3. Inspecting the outcome
# Now that we applied the GRAC shift we can observe the density and see how well we do.

# +
grac_ip = -1 * grac_data["eigenvalues"].np[dft_wfn.nalpha()] * psi4.constants.hartree2ev

println("Neon Ionization Potential (eV)")
printfmt("SCF     {:10.4f}\n", scf_ip)
printfmt("PBE0    {:10.4f}\n", dft_ip)
printfmt("PBE0-AC {:10.4f}\n", grac_ip)
printfmt("CISD    {:10.4f}\n", ci_ip)
printfmt("NIST    {:10.4f}\n", 21.5645)
# -

# We observe that our ionization potenital is now much closer to the correct ionization potenital and while the density is not perfect, it is much closer to the CISD density.

R, data = one_electron_radial([scf_wfn.Da(), dft_wfn.Da(), grac_data["Da"], ci_wfn.Da()], Vpot)
data_dict = Dict("SCF"=> data[1],
                 "DFT"=> data[2],
                 "DFT-AC"=> data[3])
plot_scatter_ratio(data_dict, data[4], R)

# Refs:
#  - M. Gruning, O. V. Gritsenko, S. J. A. van Gisbergen, and E. J. BaerendsJ. Chem. Phys., **2001**, *114*, 652


text:
export MinkowskiSumArray,
       array

# ==================================
# Minkowski sum of an array of sets
# ==================================

"""
   MinkowskiSumArray{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the Minkowski sum of a finite number of convex sets.

### Fields

- `array` -- array of convex sets

### Notes

This type assumes that the dimensions of all elements match.

The `ZeroSet` is the neutral element and the `EmptySet` is the absorbing element
for `MinkowskiSumArray`.

Constructors:

- `MinkowskiSumArray(array::Vector{<:LazySet})` -- default constructor

- `MinkowskiSumArray([n]::Int=0, [N]::Type=Float64)`
 -- constructor for an empty sum with optional size hint and numeric type
"""
struct MinkowskiSumArray{N<:Real, S<:LazySet{N}} <: LazySet{N}
   array::Vector{S}
end

isoperationtype(::Type{<:MinkowskiSumArray}) = true
isconvextype(::Type{MinkowskiSumArray{N, S}}) where {N, S} = isconvextype(S)

# constructor for an empty sum with optional size hint and numeric type
function MinkowskiSumArray(n::Int=0, N::Type=Float64)
   arr = Vector{LazySet{N}}()
   sizehint!(arr, n)
   return MinkowskiSumArray(arr)
end

# ZeroSet is the neutral element for MinkowskiSumArray
@neutral(MinkowskiSumArray, ZeroSet)

# EmptySet and Universe are the absorbing elements for MinkowskiSumArray
@absorbing(MinkowskiSumArray, EmptySet)
# @absorbing(MinkowskiSumArray, Universe)  # TODO problematic

# add functions connecting MinkowskiSum and MinkowskiSumArray
@declare_array_version(MinkowskiSum, MinkowskiSumArray)

"""
   array(msa::MinkowskiSumArray{N, S}) where {N<:Real, S<:LazySet{N}}

Return the array of a Minkowski sum of a finite number of convex sets.

### Input

- `msa` -- Minkowski sum array

### Output

The array of a Minkowski sum of a finite number of convex sets.
"""
function array(msa::MinkowskiSumArray{N, S}) where {N<:Real, S<:LazySet{N}}
   return msa.array
end

"""
   dim(msa::MinkowskiSumArray)

Return the dimension of a Minkowski sum of a finite number of sets.

### Input

- `msa` -- Minkowski sum array

### Output

The ambient dimension of the Minkowski sum of a finite number of sets.
"""
function dim(msa::MinkowskiSumArray)
   return length(msa.array) == 0 ? 0 : dim(msa.array[1])
end

"""
   σ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}

Return the support vector of a Minkowski sum of a finite number of sets in a
given direction.

### Input

- `d`   -- direction
- `msa` -- Minkowski sum array

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the summand sets.
"""
function σ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}
   return σ_helper(d, msa.array)
end

"""
   ρ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}

Return the support function of a Minkowski sum array of a finite number of sets
in a given direction.

### Input

- `d`   -- direction
- `msa` -- Minkowski sum array

### Output

The support function in the given direction.

### Algorithm

The support function of the Minkowski sum of sets is the sum of the support
functions of each set.
"""
function ρ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}
   return sum([ρ(d, Xi) for Xi in msa.array])
end

"""
	isbounded(msa::MinkowskiSumArray)

Determine whether a Minkowski sum of a finite number of convex sets is bounded.

### Input

- `msa` -- Minkowski sum of a finite number of convex sets

### Output

`true` iff all wrapped sets are bounded.
"""
function isbounded(msa::MinkowskiSumArray)
   return all(isbounded, msa.array)
end

"""
   isempty(msa::MinkowskiSumArray)

Return if a Minkowski sum array is empty or not.

### Input

- `msa` -- Minkowski sum array

### Output

`true` iff any of the wrapped sets are empty.
"""
function isempty(msa::MinkowskiSumArray)
   return any(isempty, array(msa))
end

function concretize(msa::MinkowskiSumArray)
    a = array(msa)
    @assert !isempty(a) "an empty Minkowski sum is not allowed"
    X = msa
    @inbounds for (i, Y) in enumerate(a)
        if i == 1
            X = concretize(Y)
        else
            X = minkowski_sum(X, concretize(Y))
        end
    end
    return X
end


text:
"""
    The Tower of Hanoi

# Brief:
    The Tower of Hanoi is a mathematical game or puzzle consisting of three rods and a number
    of disks of various diameters, which can slide onto any rod.
    The puzzle begins with the disks stacked on one rod in order of decreasing size,
    the smallest at the top, thus approximating a conical shape.
    The objective of the puzzle is to move the entire stack to the last rod,
    obeying the following rules:
        1. Only one disk may be moved at a time.
        2. Each move consists of taking the upper disk from one of the stacks and placing it on
           top of another stack or on an empty rod.
        3. No disk may be placed on top of a disk that is smaller than it.

# Complexity: O(2^n)

# Functions
    - solveUtil(tower1, tower2, tower3, n::Int, solution::Array{Pair}) - The recursive function which finds the solution and saves it to
                                                                         the array of pairs where each pair represents one move
    - solve(tower1, tower2, tower3,n::Int) - The main function for finding the set of steps to move n rings from tower 1 to tower 3, using tower 2
    - printSolution(solution::Array{Pair}) - Prints the solution given as an array of pairs by printing each pair ( step ) in a different line

# Reference:- [Wikipedia](https://en.wikipedia.org/wiki/Tower_of_Hanoi)

# Contributed by:- [Nikola Mircic](https://github.com/Nikola-Mircic)
"""

module Hanoi
    function solveUtil(tower1, tower2, tower3, n::Int, solution::Array{Pair})
        if n==1
            push!(solution, tower1 => tower3) # There is only one ring, so just move it from tower1 to tower3
        else
            solveUtil(tower1, tower3, tower2, n-1, solution) # Move n-1 rings from tower1 to tower2 to free the bottom ring
            push!(solution, tower1 => tower3) # Move the bottom ring from tower1 to tower3
            solveUtil(tower2, tower1, tower3, n-1, solution) # Now move those n-1 rings from tower2 to tower3
        end
    end

    function solve(tower1, tower2, tower3,n::Int)
        solution = Array{Pair}(undef, 0)
        solveUtil(tower1, tower2, tower3, n, solution)
        return solution
    end

    function printSolution(solution::Array{Pair})
        for step in solution
            println(step)
        end
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
"""
    ismeanfeasible(m::Model; [toler::Number = 1e-5, solver = Clp.ClpSolver()])
Returns a `Bool` indicating if the mean stored in `FlexibilityData.RVmeans` is feasible, meaning that it
lies inside the feasible region. This check is done using the so-called feasibility function.

**Arguments**
- `m::Model` The flexibility model.

**Keyword Arguments**
- `toler::Number = 1e-5` The numerical tolerance for checking the feasibility.
- `solver = Clp.ClpSolver()` The solver, any LP or NLP solver shoudl work.

```julia
julia> ismeanfeasible(m)
true

julia> ismeanfeasible(m, solver = GurobiSolver(OutputFlag = 0))
Academic license - for non-commercial use only
true
```
"""
function ismeanfeasible(m::Model; toler::Number = 1e-5, solver = Clp.ClpSolver())
    # Make the input dictionary
    input_dict = MakeInputDict(m)

    # Parse the necessary data dimensions
    n_z = input_dict["dims"][1]
    n_θ = input_dict["dims"][2]
    n_x = input_dict["dims"][3]
    n_f = input_dict["dims"][4]
    n_h = input_dict["dims"][5]

    # Prepare the JuMP model
    m_solve = Model(solver = solver)

    # Initialize variables
    @variable(m_solve, θ[1:n_θ])
    @variable(m_solve, u)
    if n_z != 0
        @variable(m_solve, z[1:n_z])
    end
    if n_x != 0
        @variable(m_solve, x[1:n_x])
    end

    # Set the objective
    @objective(m_solve, Min, u)

    # Fix θ and set values
    for i = 1:n_θ
        JuMP.fix(θ[i], input_dict["theta_nom"][i])
    end

    # Add the constraints
    exprs = AddSystemExpressions(m_solve, input_dict)
    @constraint(m_solve, f[j = 1:n_f], exprs[1][j] <= u)
    if n_h != 0
        @constraint(m_solve, h[i = 1:n_h], exprs[2][i] == 0)
    end

    # Solve and return
    status = solve(m_solve)
    return getvalue(u) <= toler
end

"""
    ComputeCenter(m::Model, center::Symbol, solver, toler::Number, only_positive::Bool)
Returns a center point that can be used to replace the mean if desired.

**Arguments**
- `m::Model` The flexibility model.
- `center::Symbol` Indicates the type of center, accepted arguments are `:feasible` and `:analytic`.
- `solver` The solver which must be an NLP solver for the analytic center.
- `toler::Number` The tolerance to check solution validity.
- `only_positive::Bool` Indicates if the center need by strictly positive.
"""
function ComputeCenter(m::Model, center::Symbol, solver, toler::Number, only_positive::Bool)
    # Check input
    if center != :analytic && center != :feasible
        error("Invalid argument center = $center, must be :analytic or :feasible")
    end

    # Make the input dictionary
    input_dict = MakeInputDict(m)

    # Parse the necessary data dimensions
    n_z = input_dict["dims"][1]
    n_θ = input_dict["dims"][2]
    n_x = input_dict["dims"][3]
    n_f = input_dict["dims"][4]
    n_h = input_dict["dims"][5]

    # Prepare the JuMP model
    m_solve = Model(solver = solver)

    # Initialize variables
    if only_positive
        @variable(m_solve, θ[1:n_θ] >= 0)
    else
        @variable(m_solve, θ[1:n_θ])
    end
    if center == :feasible
        @variable(m_solve, u)
    else
        @variable(m_solve, s[1:n_f] >= 0)
    end
    if n_z != 0
        @variable(m_solve, z[1:n_z])
    end
    if n_x != 0
        @variable(m_solve, x[1:n_x])
    end

    # Set objective
    if center == :feasible
        @objective(m_solve, Min, u)
    else
        @NLobjective(m_solve, Min, -sum(log(s[i]) for i = 1:n_f))
    end

    # Add constraints
    exprs = AddSystemExpressions(m_solve, input_dict)
    if center == :feasible
        @constraint(m_solve, f[j = 1:n_f], exprs[1][j] <= u)
    else
        @constraint(m_solve, f[j = 1:n_f], exprs[1][j] + s[j] == 0)
    end
    if n_h != 0
        @constraint(m_solve, h[i = 1:n_h], exprs[2][i] == 0)
    end

    # Solve and return
    status = solve(m_solve)
    if center == :feasible && getvalue(u) > toler
        upper_b = getvalue(u)
        @warn "Optimized mean is not feasible, can only achieve inequalities with upper bound u = $upper_b"
    end
    return getvalue(θ)
end

"""
    findcenteredmean(m::Model; [center::Symbol = :feasible, solver = Ipopt.IpoptSolver(print_level = 0), toler::Number = 1e-5, update_mean::Bool = false, only_positive::Bool = false])
Returns a center point based on the analytic or feasible center. The result can overwrite the mean stored in `FlexibilityData.RVmeans`
if desired. This is a wrapper function for [`ComputeCenter`](@ref).

**Arguments**
- `m::Model` The flexibility model.

**Keyword Arguments**
- `center::Symbol = :feasible` Indicates the type of center, accepted arguments are `:feasible` and `:analytic`.
- `solver = Ipopt.IpoptSolver(print_level = 0)` The solver which must be an NLP solver for the analytic center.
- `toler::Number = 1e-5` The tolerance to check solution validity.
- `update_mean::Bool = false` Indicates if the computed center should overwrite `FlexibilityData.RVmeans`.
- `only_positive::Bool = false` Indicates if the center need by strictly positive.

```julia
julia> findcenteredmean(m, only_positive = true)
4-element Array{Float64,1}:
 1684.74
   79.0718
  195.073
    0.0

julia> findcenteredmean(m, center = :analytic, update_mean = true)
4-element Array{Float64,1}:
  898.125
 -507.214
  594.544
  317.23
```
"""
function findcenteredmean(m::Model; center::Symbol = :feasible, solver = Ipopt.IpoptSolver(print_level = 0), toler::Number = 1e-5, update_mean::Bool = false, only_positive::Bool = false)
    center = ComputeCenter(m, center, solver, toler, only_positive) # included in functions.jl
    if update_mean
        flex_data = getflexibilitydata(m)
        flex_data.RVmeans = center
    end
    return center
end


text:
using Plots

eff = 0.15
area = 30

teta_p = deg2rad(170)
alfa_p = deg2rad(35)

lat = deg2rad(59.3)

delta_s(n) = deg2rad(-23.44) * cos((2pi / 365) * n)

omega(t) = deg2rad((15 * t) - 180)

alfa_s(n, t) = asin(sin(lat) * sin(delta_s(n)) + cos(lat) * cos(delta_s(n)) * cos(omega(t)))

teta_s(n, t) = pi + (omega(t) / abs(omega(t))) * acos((sin(lat) * sin(alfa_s(n, t)) - sin(delta_s(n))) / (cos(lat) * cos(alfa_s(n, t))))

I_0 = 1360
I(n, t) = 1.1 * I_0 * 0.7^((1 / (abs(sin(alfa_s(n, t)))+sin(alfa_s(n, t)))/2)^0.678)  # underlig abs medelvärde för att få 0 om negativt värde
# Detta ser till att vi får reella talet 0 när solen inte lyser på våran platts på jorden.

I_p(n, t) = I(n, t) * (cos(teta_p - teta_s(n, t)) * cos(alfa_p - alfa_s(n, t)) + sin(alfa_s(n, t)) * sin(alfa_p) * (1 - cos(teta_p - teta_s(n, t))))

weather_loss = 1-(1821/(365*12))

P(x) = (*(((abs(I_p(floor(x/24), x))+I_p(floor(x/24), x))/2),weather_loss , eff, area))

s_date = 152
e_date = 181
step = 0.1

x = 24*s_date:step:24*(e_date)+24
x = filter(x->x%12!=0,x)
# samma som i sun, men nu behövs det faktiskt eftersom vi kollar på ett större spann och behöver mer presicion.
# vi kan inte häller använda plot funktionen eftersom vi integrerar nu.

kwh = 0

for i in x
    global kwh += i*P(i)
end

print(kwh)


text:
## Static server from Meddle.jl, extended to serve index.html on a request to root

using URIParser: parse_url
using Meddle
using HttpServer
using HttpCommon: FileResponse

export ServeStatic

path_in_dir(p::String, d::String) = length(p) > length(d) && p[1:length(d)] == d

function strip_trailing_slash(p::String)
  if p[length(p)] == '/' && p != "/"
    return p[1:length(p)-1]
  end
  p
end

function ServeStatic(root::String)
  Midware() do req::MeddleRequest, res::Response
    m = match(r"^/+(.*)$", req.state[:resource])
    if m != nothing
      path = normpath(root, m.captures[1])
      # Protect against dir-escaping
      if !path_in_dir(path, root)
        return respond(req, Response(400)) # Bad Request
      end
      path = strip_trailing_slash(path)
      # Serve the requested file, if it exists
      if isfile(path)
        # res.data = readall(path)
        return respond(req, FileResponse(path))
      end
      # Serve index.html if a request is made to the root
      if (path == normpath(root))
        ixpath = path * "/index.html"
        if isfile(ixpath)
          res = FileResponse(ixpath)
          return respond(req, res)
        end
      end
      req, res
    end
  end
end


text:
using WAV
using AutoGainCtrl

wav, fs = wavread("./doc/123.wav")

new = agc(wav,fs,maxvalue=0.5,minstep=-0.99); # deepcopy wav and then operate in-place
agc!(wav,fs,maxvalue=0.15,minstep=-0.5)       # in-place operation

wavwrite(wav, "./doc/agced1.wav", Fs=fs, nbits=32)
wavwrite(new, "./doc/agced2.wav", Fs=fs, nbits=32)

max_stars_repo_name: sonosole/AGC

max_issues_repo_name: sonosole/AutoGainCtrl

max_forks_repo_name: sonosole/AutoGainCtrl


text:
using DashBootstrapComponents, DashHtmlComponents

card = dbc_card(
    dbc_cardbody([
        html_h4("Title", className = "card-title"),
        html_h6("Card subtitle", className = "card-subtitle"),
        html_p(
            "Some quick example text to build on the card title and make " *
            "up the bulk of the card's content.",
            className = "card-text",
        ),
        dbc_cardlink("Card link", href = "#"),
        dbc_cardlink("External link", href = "https://google.com"),
    ]),
    style = Dict("width" => "18rem"),
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export Arc, get_ilabel, get_olabel, get_weight, get_nextstate
"""
  `Arc(ilab::Int,olab::Int,w::W,n::Int)`

Constructs an arc with input label `ilab`, output label `olab`, weight 'weight' and nextstate `n`.
`ilab` and `olab` must be integers consistent with a input/output symbol table of the WFST at use.

Mainly used for internals see [`add_arc!`](@ref) for a simpler way of adding arcs to a WFST.
"""
struct Arc{W,D}
  ilabel::D      # input label
  olabel::D      # output label
  weight::W      # weight
  nextstate::D   # destination id
end

import Base: ==

function ==(a::Arc,b::Arc; check_weight=true)
  i = (a.ilabel == b.ilabel)
  o = (a.olabel == b.olabel)
  d = (a.nextstate == b.nextstate)
  if check_weight
    w = (a.weight == b.weight)
    return i && o && d && w
  else
    return i && o && d
  end
end

function Base.show(io::IO, A::Arc, isym::Dict, osym::Dict)
  ilabel = idx2string(A.ilabel,isym)
  olabel = idx2string(A.olabel,osym)
  print(io, "$(ilabel):$(olabel)/$(A.weight) → ($(A.nextstate))")
end

function Base.show(io::IO, A::Arc)
  print(io, "$(A.ilabel):$(A.olabel)/$(A.weight) → ($(A.nextstate))")
end

function idx2string(label,sym)
  if iseps(label)
    return "ϵ"
  else
    return "$(sym[label])"
  end
end

"""
`get_ilabel(A::Arc)`

Returns the input label index of the arc `A`.
"""
get_ilabel(A::Arc) = A.ilabel

"""
`get_olabel(A::Arc)`

Returns the input label index of the arc `A`.
"""
get_olabel(A::Arc) = A.olabel

"""
`get_weight(A::Arc)`

Returns the weight the arc `A`.
"""
get_weight(A::Arc) = A.weight

"""
`get_nextstate(A::Arc)`

Returns the state that the arc `A` is pointing to.
"""
get_nextstate(A::Arc) = A.nextstate

change_next(a::Arc{W,D}, nextstate::D) where {W,D} = Arc{W,D}(get_ilabel(a),
                                                              get_olabel(a),
                                                              get_weight(a),
                                                              nextstate)
Base.inv(a::Arc) = Arc(get_olabel(a),get_ilabel(a),get_weight(a),get_nextstate(a))
proj(a::Arc, get_iolabel::Function) = Arc(get_iolabel(a),get_iolabel(a),get_weight(a),get_nextstate(a))


text:
# States for time correlated noise

type q_noise
	q_SW::Float64
	q_LW::Float64
	q_P::Float64
	q_Ta::Float64
	q_RH::Float64
	q_Ua::Float64
end

# Perturb input data

function perturb_input(q_vec, ipart, timestep, SW, LW, P, Ta, RH, Ua)

  # Shortwave radiation - random errors

  decorr_SW = 3.0

  error_SW = min(SW, 109.1)

  q_vec[ipart].q_SW = timecorr_noise(q_vec[ipart].q_SW, timestep, decorr_SW)

  SW_noise = error_SW * q_vec[ipart].q_SW

  SW = SW + SW_noise

  SW = max(0.0, SW)

  # Longwave radiation - random errors

  decorr_LW = 4.7

  error_LW = 20.8

  q_vec[ipart].q_LW = timecorr_noise(q_vec[ipart].q_LW, timestep, decorr_LW)

  LW_noise = error_LW * q_vec[ipart].q_LW

  LW = LW + LW_noise

  # Precipitation - random errors

  decorr_P = 2.0

  mu_P = -0.19
  sigma_P = 0.61

  q_vec[ipart].q_P = timecorr_noise(q_vec[ipart].q_P, timestep, decorr_P)

  P_noise = exp(mu_P + sigma_P * q_vec[ipart].q_P)

  P = max(P * P_noise, 0.0)

  # Air temperature - random errors

  decorr_Ta = 4.8

  error_Ta = 0.9

  q_vec[ipart].q_Ta = timecorr_noise(q_vec[ipart].q_Ta, timestep, decorr_Ta)

  Ta_noise = error_Ta * q_vec[ipart].q_Ta

  Ta = Ta + Ta_noise

  # Relative humidity - random noise

  decorr_RH = 8.4

  error_RH = 8.9

  q_vec[ipart].q_RH = timecorr_noise(q_vec[ipart].q_RH, timestep, decorr_RH)

  RH_noise = error_RH * q_vec[ipart].q_RH

  RH = RH + RH_noise

  RH = max(0.0, RH)
  RH = min(100.0, RH)

  # Wind speed - random noise

  decorr_Ua = 8.2

  mu_Ua = -0.14
  sigma_Ua = 0.53

  q_vec[ipart].q_Ua = timecorr_noise(q_vec[ipart].q_Ua, timestep, decorr_Ua)

  Ua_noise = exp(mu_Ua + sigma_Ua * q_vec[ipart].q_Ua)

  Ua = Ua * Ua_noise

  Ua = max(0.5, Ua)
  Ua = min(25.0, Ua)

  return SW, LW, P, Ta, RH, Ua

end


text:
# Define variables used throughout the package and functions to initialize them.

export GainProfile
export gen_gain_2lv, hole_burning!

## Two-level atoms
# Generate the default gain curve of SALT, which describes two-level atoms.
# For multiple classes of atoms describing inhomogeneous broadening, for which ωₐ and γ˔ are
# vectors, use gen_gain_2lv.(ωₐ, γ˔).
gen_gain_2lv(ωₐ::Real, γ˔::Real) = ω::Number -> γ˔ / (ω - ωₐ + im * γ˔)  # scalar
gen_gain′_2lv(ωₐ::Real, γ˔::Real) = ω::Real -> -γ˔ / (ω - ωₐ + im * γ˔)^2  # scalar

# Below, we define the formula for |γ|².  This can be easily calculated from γ, but the
# users may accidentally evaluate γ(ω) for complex ω and take its squared absolute value.
# That will give a wrong result, because the imaginary part of ω will be added to iγ in the
# denominator when calculating the absolute value.  To prevent this unfortunate accident
# from happening, we define a formula for |γ|² that accepts only real ω and ask the users to
# use it.
gen_abs2gain_2lv(ωₐ::Real, γ˔::Real) = ω::Real -> γ˔^2 / ((ω - ωₐ)^2 + γ˔^2)  # scalar; note ω is real for this
gen_abs2gain′_2lv(ωₐ::Real, γ˔::Real) = ω::Real -> -2γ˔^2 * (ω - ωₐ) / ((ω - ωₐ)^2 + γ˔^2)^2  # scalar; note ω is real for this

# Below, allow vectors and matrices to be PETSc ones if their sizes are 3×(# of grid points).
# If their sizes are the number of modes, keep them Julia vectors.

# About avoiding overtyping container-type variables:
#
# Below, for read-only container-type variables, I specify neither the container type nor
# eltype concretely.  [What are the examples of read-only container-type variables?]
#
# However, for variables into which I write, again I don't specify the container type
# concretely (so that I can allow various container type, such as Julia vectors and PETSc
# vectors), but I do specify the eltype, because otherwise I get InexactError that complains
# about incompatible element type (e.g., when attempting to store ComplexF64 in Float64).
#
# The same principle applies to writing functions as well as types.  If the function changes
# the contents of some variable, that variable is where I write, so its eltype should be
# specified to avoid InexactError.  For example, in create_A!, the argument A is where the
# output is stored, so it had better have ComplexF64 as the eltype, so I type it as
# AbsMatComplexF.  On the other hand, the argument CC is read-only, so it is typed AbsMatNumber.

# Parameters defining the SALT problem
# Consider including CC to gp, if I am really going to use ωₐ for PML for all modes.
mutable struct GainProfile{VF<:AbsVecFloat,G,G′,AG,AG′,W}  # VF can be PETSc vectors
    gain::Vector{G}  # gain curve
    gain′::Vector{G′}  # derivative of gain curve
    abs2gain::Vector{AG}  # squared absolute value of gain curve
    abs2gain′::Vector{AG′}  # derivative of squared absolute value of gain curve
    D₀::VF  # pump strength; length(D₀) = number of DOFs
    wt::Vector{W}  # function that calculates contribution of D₀ to kth atomic class; D₀ₖ = wt.(D₀)

    function GainProfile{VF,G,G′,AG,AG′,W}(gain::AbsVec{G},
                                           gain′::AbsVec{G′},
                                           abs2gain::AbsVec{AG},
                                           abs2gain′::AbsVec{AG′},
                                           D₀::AbsVecReal,
                                           wt::AbsVec{W}) where {VF<:AbsVecFloat,G,G′,AG,AG′,W}
        K = length(wt)  # number of atomic classes
        length(gain)==length(gain′)==length(abs2gain)==length(abs2gain′)==K ||
            @error "length.((gain,gain′,abs2gain,abs2gain′)) = $(length.((gain,gain′,abs2gain,abs2gain′))) should be all same."

        return new(gain, gain′, abs2gain, abs2gain′, D₀, wt)
    end
end

function GainProfile(gain::AbsVec, gain′::AbsVec, abs2gain::AbsVec, abs2gain′::AbsVec, D₀::AbsVecReal,
                     wt::AbsVec=(K=length(gain); [(d::Real->d/K) for k=1:K]))  # default wt: even distribution
    D₀_new = similar(D₀,Float)
    copyto!(D₀_new, D₀)

    return GainProfile{typeof(D₀_new),eltype.((gain,gain′,abs2gain,abs2gain′,wt))...}(gain, gain′, abs2gain, abs2gain′, D₀_new, wt)
end

# Convenience constructors for single atomic class
GainProfile(gain::Any, gain′::Any, abs2gain::Any, abs2gain′::Any, D₀::AbsVecReal) =
    GainProfile([gain], [gain′], [abs2gain], [abs2gain′], D₀)
GainProfile(gain::Any, gain′::Any, abs2gain::Any, abs2gain′::Any, N::Integer) =
    GainProfile(gain, gain′, abs2gain, abs2gain′, VecFloat(undef,N))

# Constructors using automatic differentiation
function GainProfile(gain::AbsVec,
                     D₀::AbsVecReal,
                     wt::AbsVec=(K=length(gain); [(d::Real->d/K) for k=1:K]))  # default wt: even distribution
    gain′ = [ω->(g(x) = derivative(x->real(γ(x+ω)),x) + im * derivative(x->imag(γ(x+ω)),x); g(0)) for γ = gain]
    abs2gain = [ω->abs2(γ(ω)) for γ = gain]
    abs2gain′ = [ω->derivative(abs2γ, ω) for abs2γ = abs2gain]

    return GainProfile(gain, gain′, abs2gain, abs2gain′, D₀, wt)
end

# Convenience constructors for single atomic class
GainProfile(gain::Any, D₀::AbsVecReal) = GainProfile([gain], D₀)
GainProfile(gain::Any, N::Integer) = GainProfile(gain, VecFloat(undef,N))

Base.length(gp::GainProfile) = length(gp.wt)  # number of atomic classes

# Evaluate 1 + hole-burning term = 1 + ∑|γaψ|².
function hole_burning!(hb::AbsVecNumber,  # output
                       abs2gain::Any,  # function that takes ω and produces |γ(ω)|²
                       ω::AbsVecReal,  # vector of lasing frequencies
                       a²::AbsVecReal,  # vector of squared amplitudes of unnormalized eigenmodes
                       abs2ψ::AbsVec{<:AbsVecReal})  # vector of absolute squares of normalized eigenmodes
    hb .= 1  # initialize
    for m = 1:length(a²)
        if a²[m] ≠ 0
            # @info "a²[$m] = $(a²[m]), ‖ψ[$m]‖ = $(norm(ψ[m]))"
            hb .+=  (abs2gain(ω[m]) * a²[m]) .* abs2ψ[m]
        end
    end
    # @info "‖hb‖ = $(norm(hb))"

    return nothing
end


text:
using StanSample, Distributions

N = 100;
df = DataFrame(
    height = rand(Normal(10, 2), N),
    leg_prop = rand(Uniform(0.4, 0.5), N),
);
df.leg_left = df.leg_prop .* df.height + rand(Normal(0, 0.02), N);
df.leg_right = df.leg_prop .* df.height + rand(Normal(0, 0.03), N);

stan6_1 = "
data {
  int <lower=1> N;
  vector[N] H;
  vector[N] LL;
  vector[N] LR;
}
parameters {
  real a;
  vector[2] b;
  real <lower=0> sigma;
}
model {
  vector[N] mu;
  mu = a + b[1] * LL + b[2] * LR;
  a ~ normal(10, 100);
  b ~ normal(2, 10);
  sigma ~ exponential(1);
  H ~ normal(mu, sigma);
}
";

m6_1s = SampleModel("m6.1s", stan6_1);
data = (H = df.height, LL = df.leg_left, LR = df.leg_right, N = size(df, 1));

rc6_1s = stan_sample(m6_1s; data);
if success(rc6_1s)
    nt6_1s = read_samples(m6_1s)
end

nt6_1s.b |> display
println()
mean(nt6_1s.b, dims=2) |> display

init = (a = 2.0, b = [1.0, 2.0], sigma = 1.0)
rc6_2s = stan_sample(m6_1s; data, init);
if success(rc6_1s)
    nt6_2s = read_samples(m6_1s)
end

nt6_2s.b |> display
println()
mean(nt6_2s.b, dims=2) |> display
println()

read_summary(m6_1s) |> display
println()

post6_1s_df = read_samples(m6_1s; output_format=:dataframe) |> display
println()

part6_1s = read_samples(m6_1s; output_format=:particles) |> display


text:
__precompile__(true)

"""
The module for [Mersenne Twisters](@ref).

Currently only provide one RNG type:

- [`MT19937`](@ref)
"""
module MersenneTwisters

export MT19937

include("bases.jl")
include("main.jl")

end


text:
#
# generate Artificial Data on several manifolds
#
#
#
import LinearAlgebra: I, Diagonal
export artificialS1Signal, artificialS1SlopeSignal, artificialInSARImage
export artificialSPDImage, artificialSPDImage2
export artificialS2WhirlImage, artificialS2WhirlPatch
export artificialS2RotationsImage
export artificialS2WhirlPatch, artificialS2Lemniscate

@doc doc"""
    artificialInSARImage([pts=500, pointType=S1Point])
generate an artificial InSAR image, i.e. phase valued data, of size `pts` x
`pts` points.
The `pointType` – (`S1Point`) provide the point type, i.e. manifold the data lives
  on. Possible values: [`S1Point`](@ref), [`SnPoint`](@ref), [`RnPoint`](@ref),
  where for the second type, the data is embeded (point wise) in $\mathbb R^2$.
"""
function artificialInSARImage(pts::Integer,
      pointType::Union{Type{S1Point},Type{RnPoint},Type{SnPoint},Type{Real}}=S1Point)
  # variables
  # rotation of ellipse
  aEll = 35.0; cosE = cosd(aEll); sinE = sind(aEll)
  aStep = 45.0; cosA = cosd(aStep); sinA = sind(aStep)
  # main and minor axis of the ellipse
  axes_inv = [6,25]
  # values for the hyperboloid
  midPoint = [0.275;0.275]
  radius = 0.18
  values = [range(-0.5,0.5,length = pts)...]
  # Steps
  aSteps = 60.0; cosS = cosd(aSteps); sinS = sind(aSteps)
  l = 0.075
  midP = [-.475,-.0625];#.125, .55]
  img = zeros(Float64,pts,pts)
  for j in eachindex(values), i in eachindex(values)
    # ellipse
    Xr = cosE*values[i] - sinE*values[j]
    Yr = cosE*values[j] + sinE*values[i]
    v = axes_inv[1]*Xr^2 + axes_inv[2]*Yr^2
    k1 = v <= 1.0 ? 10.0*pi*Yr : 0.0
    # circle
    Xr = cosA*values[i] - sinA*values[j]
    Yr = cosA*values[j] + sinA*values[i]
    v = ( (Xr-midPoint[1])^2 + (Yr-midPoint[2])^2 )/radius^2
    k2 = v <= 1.0 ? 4.0*pi*(1.0-v) : 0.0
    #
    Xr = cosS*values[i] - sinS*values[j]
    Yr = cosS*values[j] + sinS*values[i]
    k3=0.0
    for m=1:8
      k3 += (abs(Xr+midP[1]+m*l)+abs(Yr+midP[2]+m*l)) <= l ? 2*pi*(m/8) : 0.0
    end
    img[i,j] = mod(k1+k2+k3-pi,2*pi)+pi
  end
  if isa(pointType,Type{SnPoint}) # SnPoint: wrap & embed
    return embed.(Ref(Circle()),S1Point.(symRem.(img)))
  elseif isa(pointType,Type{S1Point}) # S1Point: wrap
    return S1Point.(symRem.(img))
  else
    return pointType.(img)
  end
end

@doc doc"""
    artificialS1SlopeSignal([pts=500, slope=4.])

Creates a Signal of (phase-valued) data represented on the
[`Circle`](@ref)` `[`Manifold`](@ref) with increasing slope.

# Optional
* `pts` – (`500`) number of points to sample the function.
* `slope` – (`4.0`) initial slope that gets increased afterwards
"""
function artificialS1SlopeSignal(pts::Integer = 500, slope::Float64=4.)
    t = range(0., 1., length=pts)
    f = zero(t)
    f[ t .<= 1/6] .= -π/2 .+ slope * π/8 * t[ t .<= 1/6]
    # In the following terms, the first max
    f[(1/6 .< t) .& (t .<= 1/3)] .= max(f[f .!= 0]...) .- slope * π/4 * 1/6 .+ slope * π/4 .* t[ (1/6 .< t) .& ( t.<=1/3 )]
    f[(1/3 .< t) .& (t .<= 1/2)] .= max(f[f .!= 0]...) .- slope * π/2 * 1/3 .+ slope * π/2 * t[ (1/3 .< t) .& ( t .<= 1/2)]
    f[(1/2 .< t) .& (t .<= 2/3)] .= max(f[f .!= 0]...) .- slope * π * 1/2 .+ slope * π * t[(1/2 .< t) .& (t .<= 2/3)]
    f[(2/3 .< t) .& (t .<= 5/6)] .= max(f[f .!= 0]...) .- slope * 2*π * 2/3 .+ slope * 2 * π * t[(2/3 .< t) .& (t .<= 5/6)]
    f[ 5/6 .< t] .= max(f[f .!= 0]...) .- slope * 4 * π * 5/6 .+ slope * 4 * π * t[ 5/6 .< t]
    return S1Point.(symRem.(f))
end

@doc doc"""
    artificialS1Signal([pts=500,pointType=S1Point])

generate a real-valued signal having piecewise constant, linear and quadratic
intervals with jumps in between. If the resulting manifold the data lives on,
is the [`Circle`](@ref) (the default [`S1Point`](@ref)) the data is also wrapped
to $[-\pi,\pi)$.

# Optional
* `pts` – (`500`) number of points to sample the function
* `pointType` – (`S1Point`) provide the point type, i.e. manifold the data lives
  on. Possible values: [`S1Point`](@ref), [`SnPoint`](@ref), [`RnPoint`](@ref),
  where for the second type, the data is embeypeded (point wise) in $\mathbb R^2$.
"""
function artificialS1Signal(pts::Integer=500,
    pointType::Union{Type{S1Point},Type{RnPoint},Type{SnPoint},Type{Real}}=S1Point)
  t = range(0., 1., length=pts)
  f = artificialS1Signal.(t)
  if isa(pointType,Type{SnPoint}) # SnPoint: wrap & embed
    return embed.(Ref(Circle()),S1Point.(symRem.(f)))
  elseif isa(pointType,Type{S1Point}) # S1Point: wrap
    return S1Point.(symRem.(f))
  else
    return pointType.(f)
  end
end
@doc doc"""
    artificialS1Signal(x)
evaluate the example signal $f(x), x\in [0,1]$,
of phase-valued data introduces in Sec. 5.1 of

> Bergmann, Laus, Steidl, Weinmann, Second Order Differences of Cyclic Data and
> Applications in Variational Denoising, SIAM J. Imaging Sci., 7(4), 2916–2953, 2014.
> doi: [10.1137/140969993](https://dx.doi.org/10.1137/140969993)

for values outside that intervall, this Signal is `missing`.
"""
function artificialS1Signal(x::Number)
    if x < 0
        y = missing
    elseif x <= 1/4
        y = - 24 * π * (x-1/4)^2  +  3/4 * π
    elseif x <= 3/8
        y = 4 * π * x  -  π / 4
    elseif x <= 1/2
        y = - π * x - 3*π/8;
    elseif x <= (3*0 + 19)/32
        y = - (0+7)/8*π
    elseif x <= (3*1 + 19)/32
        y = - (1+7)/8*π
    elseif x <= (3*2 + 19)/32
        y = - (2+7)/8*π
    elseif x <= (3*3 + 19)/32
        y = - (3+7)/8*π
    elseif x <= 1
        y = 3 / 2 * π * exp(8-1/(1-x))  -  3/4*π
    else
        y = missing
    end
    return y
end
@doc doc"""
    artificialS2WhirlImage([pts=64])
generate an artificial image of data on the 2 sphere,

# Arguments
* `pts` – (`64`) size of the image in `pts`$\times$`pts` pixel.
"""
function artificialS2WhirlImage(pts::Int=64)
  M = Sphere(2)
  N = Power(M, (pts,pts) )
  # background - default rotations
  img = artificialS2RotationsImage(pts, (0.5,0.5) )
  # Set WhirlPatches
  sc = pts/64
  patchSizes = floor.( sc.* [9,9,9,9,11,11,11,15,15,15,17,21] )
  patchCenters = Integer.( floor.(
    sc.* [ [35,7] [25,41] [32,25] [7,60] [10,5] [41,58] [11,41] [23,56] [38,45] [16,28] [55,42] [51,16] ]
  ) )
  patchSigns = [1,1,-1,1,-1,-1,1,1,-1,-1,1,-1]
  for i=1:length(patchSizes)
    pS = Integer(patchSizes[i])
    pSH = Integer( floor( (patchSizes[i]-1)/2 ) )
    pC = patchCenters[:,i]
    pC = max.(Ref(1), pC.-pS) .+pS
    pSgn = patchSigns[i]
    s = pS%2==0 ? 1 : 0
    r = [ pC[1] .+ ( -pSH:(pSH+s) ), pC[2] .+ ( -pSH:(pSH+s) ) ]
    patch = artificialS2WhirlPatch(pS)
    if pSgn==-1 # opposite ?
      patch = PowPoint( opposite.(Ref(M), getValue(patch)) )
    end
    img[ r... ] = patch
  end
  return img
end
@doc doc"""
    artificialS2Rot([pts=64, rotations=(.5,.5)])
creates an image with a rotation on each axis as a parametrization.

# Optional Parameters
* `pts` – (`64`) number of pixels along one dimension
* `rotations` – (`(.5,.5)`) number of total rotations performed on the axes.
"""
function artificialS2RotationsImage(pts::Int=64,rotations::Tuple{Float64,Float64}=(.5,.5))
  M = Sphere(2)
  N = Power(M, (pts,pts) )
  img = PowPoint(Matrix{SnPoint{Float64}}(undef,pts,pts))
  north = SnPoint([1.,0.,0.])
  Rxy(a) = [cos(a) -sin(a) 0.; sin(a) cos(a) 0.; 0. 0. 1]
  Rxz(a) = [cos(a) 0. -sin(a); 0. 1. 0.; sin(a) 0. cos(a)]
  for i=1:pts
    for j=1:pts
      x = i/pts*2π*rotations[1];
      y = j/pts*2π*rotations[2];
      img[i,j] = SnPoint( Rxy(x+y)*Rxz(x-y)*[0,0,1] )
    end
  end
  return img
end

@doc doc"""
    artificialS2WhirlPatch([pts=5])

create a whirl within the `pts`$\times$`pts` patch of
[`Sphere`](@ref)`(2)`-valued image data.

# Optional Parameters
* `pts` – (`5`) size of the patch. If the number is odd, the center is the north
  pole.
"""
function artificialS2WhirlPatch(pts::Int=5)
  M = Sphere(2)
  N = Power(M, (pts,pts) )
  patch = PowPoint( fill(SnPoint([0.,0.,-1.]),pts,pts))
  scaleFactor = sqrt( (pts-1)^2 / 2 )*3/π;
  for i=1:pts
    for j=1:pts
      if i!=(pts+1)/2 || j != (pts+1)/2
        # direction within the patch
        α = atan( (j-(pts+1)/2), (i-(pts+1)/2) );
        # scaled length
        β = sqrt( (j-(pts+1)/2)^2 + (i-(pts+1)/2)^2 )/scaleFactor;
        patch[i,j] = SnPoint([sin(α)*sin(π/2-β), -cos(α)*sin(π/2-β), cos(π/2-β)])
      end
    end
  end
  return patch
end
@doc doc"""
    artificialSPDImage([pts=64, stepsize=1.5])

create an artificial image of symmetric positive definite matrices of size
`pts`$\times$`pts` pixel with a jump of size `stepsize`.
"""
function artificialSPDImage(pts::Int=64, stepsize = 1.5)
  r = range(0, stop = 1-1/pts, length=pts)
  v1 = abs.(2*pi .* r .- pi)
  v2 = pi .* r;
  v3 = range(0, stop = 3*(1-1/pts), length = 2*pts);
  data = PowPoint(fill( SPDPoint( Matrix{Float64}(I,3,3) ), pts, pts ))
  for row = 1:pts
    for col = 1:pts
      A = [cos(v1[col]) -sin(v1[col]) 0. ; sin(v1[col]) cos(v1[col]) 0. ; 0. 0. 1.]
      B = [1. 0. 0. ; 0. cos(v2[row]) -sin(v2[row]) ; 0. sin(v2[row]) cos(v2[row]) ]
      C = [ cos(v1[mod(col-row,pts)+1]) 0 -sin(v1[mod(col-row,pts)+1]);
           0. 1. 0.;
           sin(v1[mod(col-row,pts)+1]) 0. cos(v1[mod(col-row,pts)+1]) ]
      scale = [ 1 + stepsize/2 * ( (row + col) > pts ? 1 : 0)
                1 + v3[row + col] - stepsize * ( col > pts/2 ? 1 : 0)
                4 - v3[row + col] + stepsize * ( row > pts/2 ? 1 : 0) ]
      data[row, col] = SPDPoint( A * B * C * Diagonal(scale) * C' * B' * A' )
    end
  end
  return data
end
@doc doc"""
    artificialSPDImage2([pts=64, fraction=.66])

create an artificial image of symmetric positive definite matrices of size
`pts`$\times$`pts` pixel with right hand side `fraction` is moved upwards.
"""
function artificialSPDImage2(pts=64, fraction = 0.66)
  Zl = SPDPoint( 4. * Matrix{Float64}(I,3,3) )
  # create a first matrix
  α = 2. * π/3;
  β = π/3;
  B = [  1. 0. 0. ; 0. cos(β) -sin(β) ; 0. sin(β) cos(β)  ]
  A = [  cos(α) -sin(α) 0. ; sin(α) cos(α) 0. ; 0. 0. 1.  ]
  Zo = SPDPoint( A * B * Diagonal( [2., 4., 8.] ) * B' * A' )
  # create a second matrix
  α = -4. * π/3;
  β = -π/3;
  B = [  1. 0. 0. ; 0. cos(β) -sin(β) ; 0. sin(β) cos(β)  ]
  A = [  cos(α) -sin(α) 0. ; sin(α) cos(α) 0. ; 0. 0. 1.  ]
  Zt = SPDPoint( A * B * Diagonal( [ 8. / sqrt(2.), 8., sqrt(2.) ] ) * B' * A' )
  data = PowPoint(fill( SPDPoint( Matrix{Float64}(I,3,3) ), pts, pts ))
  M = SymmetricPositiveDefinite(3)
  for row = 1 : pts
    for col = 1 : pts
      # (a) from Zo a part to Zt
      C = Zo;
      if ( row > 1) # in X direction
        C = exp(M, C, log(M,C,Zt),
              (row - 1)/(2*(pts-1)) + ( ( row > fraction * pts) ? 1/2 : 0. )
            )
      end
      if ( col > 1) # and then in Y direction
        C = exp(M, C, parallelTransport(M,Zo,C, log(M,Zo,Zl)), (col - 1.)/(pts-1) )
      end
      data[row,col] = C
    end
  end
  return data
end

@doc doc"""
    artificialS2Lemniscate(p [,pts=128,a=π/2,interval=[0,2π])

generate a Signal on the [`Sphere`](@ref)`(2)` $\mathbb S^2$ by creating the
[Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli)
in the tangent space of `p` sampled at `pts` points and use `exp` to get a
signal on the [`Sphere`](@ref)`(2)`.

# Input
* `p` – the tangent space the Lemniscate is created in
* `pts` – (`128`) number of points to sample the Lemniscate
* `a` – (`π/2`) defines a half axis of the Lemniscate to cover a
   half sphere.
* `interval` – (`[0,2*π]`) range to sample the lemniscate at, the default value
  refers to one closed curve
"""
function artificialS2Lemniscate(p::SnPoint,pts::Integer=128, a::Float64=π/2.,
    interval::Array{Float64,1}=[0.,2.0*π])
    return artificialS2Lemniscate.(Ref(p),range(interval[1],interval[2],length=pts), a)
end
@doc doc"""
    artificialS2Lemniscate(p,t; a=π/2)

generate a point from the signal on the [`Sphere`](@ref)`(2)` $\mathbb S^2$ by
creating the [Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli)
in the tangent space of `p` sampled at `t` and use èxp` to obtain a point on
the [`Sphere`](@ref)`(2)`.

# Input
* `p` – the tangent space the Lemniscate is created in
* `t` – value to sample the Lemniscate at

# Optional Values
 * `a` – (`π/2`) defines a half axis of the Lemniscate to cover a
   half sphere.
"""
function artificialS2Lemniscate(p::SnPoint,t::Float64, a::Float64=π/2.)
    M = Sphere(2)
    tP = 2.0*Float64(getValue(p)[1]>=0.)-1. # Take north or south pole
    base = SnPoint([0.,0.,tP])
    xc = a * (cos(t)/(sin(t)^2+1.))
    yc = a * (cos(t)*sin(t)/(sin(t)^2+1.))
    tV = parallelTransport(M,base,p,SnTVector([xc,yc,0.]))
    return exp(M,p,tV)
end


text:
using QuantumLattices
using ExactDiagonalization.FEDW: NCBasis,bkron,matrix,boplus
using LinearAlgebra: eigen,diagm
using ExactDiagonalization.Perturbation: GSpinGenerator,jval,hamiltonianeff,projstate

lattice=Lattice("Honeycomb",[Point(PID("A",1),(0.0,0.0),(0.0,0.0)),Point(PID("B",2),[1/2,1/(2*sqrt(3))],[0.0,0.0])],vectors=[[1.0,0.0],[1/2.0,sqrt(3)/2]],neighbors=1)
config=Config{Fock{:f}}(pid->Fock{:f}(atom=pid.site, norbital=3, nspin=2, nnambu=2), lattice.pids)
table=Table(config,usualfockindextotuple)
bond = Bonds(lattice)
tij = [58.7 113.9 -7.0;113.9 58.7 -7.0; -7.0 -7.0 -194.1]

function fcmatrixob(tij::Matrix{<:Number};kwargs...)
fc = []
n,m=size(tij)

for i=1:m
    for j=1:n
        push!(fc, FockCoupling{2}(tij[j,i];orbitals=(j,i),nambus=(2,1),kwargs...) )

    end
end
res = sum(+,fc)
return res
end

tfc = fcmatrixob(tij;atoms=(2,1))
macro Lˣ_str(::String) fc"1.0im ob[2 3]" - fc"1.0im ob[3 2]" end
macro Lʸ_str(::String) fc"1.0im ob[3 1]" - fc"1.0im ob[1 3]" end
macro Lᶻ_str(::String) fc"1.0im ob[1 2]" - fc"1.0im ob[2 1]" end
macro L⁰_str(::String) fc"1.0 ob[1 1]" + fc"1.0 ob[2 2]" + fc"1.0 ob[3 3]" end
macro soc_str(::String) 0.5*(Lˣ""⊗σˣ"sp" + Lʸ""⊗σʸ"sp"+ Lᶻ""⊗σᶻ"sp") end
macro jx_str(::String) 0.5*(L⁰""⊗σˣ"sp") - Lˣ""⊗σ⁰"sp" end
macro jy_str(::String) 0.5*(L⁰""⊗σʸ"sp") - Lʸ""⊗σ⁰"sp" end
macro jz_str(::String) 0.5*(L⁰""⊗σᶻ"sp") - Lᶻ""⊗σ⁰"sp" end
macro j0_str(::String) L⁰""⊗σ⁰"sp" end

Lx = [0 0 0; 0  0 im;0 -im 0]
Ly = [0 0 -im; 0 0 0; im 0 0]
Lz = [0 im 0; -im 0 0; 0 0 0]
L0 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]

UU=4000.0
Jₕ=0.2*UU
lambda = 140.0
 t=Hopping(:t,1.0,1,couplings=tfc,amplitude=bond->((bond|>rcoord|>azimuthd ≈ 270) ? 1 : 0) )
U=Hubbard(:U,UU)
U′= InterOrbitalInterSpin(Symbol("U′"), UU-2*Jₕ)
UmJ= InterOrbitalIntraSpin(Symbol("U′-J"), UU-3*Jₕ)
J=SpinFlip(Symbol("J"), Jₕ)
Jp=PairHopping(:Jp, Jₕ)
λ=Onsite(:λ,1.0,couplings=lambda*soc"")

h0=Generator((λ,U,U′,UmJ,J,Jp),Bonds(lattice),config,half=false,table=table)
h1 =Generator((t,),Bonds(lattice),config,half=false,table=table)

h0opts = expand(h0)
h1opts = expand(h1)

basisP0=NCBasis{UInt}([6],[[(5,)]])
basisQ2=NCBasis{UInt}([6],[[(6,)]])
basisQ0=NCBasis{UInt}([6],[[(4,)]])
Pₗ0 = projstate(h0opts,basisP0,h0.table,pick=1:2)
Qₕ2 = projstate(h0opts,basisQ2,h0.table)
Qₕ0 = projstate(h0opts,basisQ0,h0.table)
Pₗij = bkron( Pₗ0, Pₗ0 )
Qₕij = boplus(bkron(Qₕ2,Qₕ0),bkron(Qₕ0,Qₕ2))
heffterm = hamiltonianeff(Pₗij, Qₕij,h1opts,h1.table)

namestri=["σᵢ⁰","σᵢˣ","σᵢʸ","σᵢᶻ"]
namsi = tuple([Symbol(nm) for nm in namestri]...)

sx=Onsite(:sx,1.0+0.0im,couplings = jx"" )
sy=Onsite(:sy,1.0+0.0im,couplings = jy"" )
sz=Onsite(:sz,1.0+0.0im,couplings = jz"" )

σx=Pₗ0.vecs'*Array(matrix(expand(sx,bond,config,false,table=table),basisP0,table))*Pₗ0.vecs
σy=Pₗ0.vecs'*Array(matrix(expand(sy,bond,config,false,table=table),basisP0,table))*Pₗ0.vecs
σz=Pₗ0.vecs'*Array(matrix(expand(sz,bond,config,false,table=table),basisP0,table))*Pₗ0.vecs
σ0=[1.0 0.0;0.0 1.0]

gsgi=GSpinGenerator(namsi,(σ0,σx,σy,σz))
namestrj=["σⱼ⁰","σⱼˣ","σⱼʸ","σⱼᶻ"]
namsj = tuple([Symbol(nm) for nm in namestrj]...)
gsgj = GSpinGenerator(namsj,(σ0,σx,σy,σz))
gsgij=kron(gsgi,gsgj)

jex2=jval(heffterm,gsgij)

txx,tyx,tzx,txy,tyy,tzy,txz,tyz,tzz = tij[1,1],tij[2,1],tij[3,1],tij[1,2],tij[2,2],tij[3,2],tij[1,3],tij[2,3],tij[3,3]
lambda1=lambda
A = -1/3*(Jₕ+9*lambda1+3*UU)/(6*Jₕ^2-UU*(3*lambda1+UU)+Jₕ*(4*lambda1+UU))
η = Jₕ/(6*Jₕ^2+(3*lambda1+UU)*(3*lambda1+2*UU)-Jₕ*(17*lambda1+8*UU))
B = 4/3*(3*Jₕ-3*lambda1-UU)/(6*Jₕ-3*lambda1-2*UU)*η

J23 = 8*A/9*(-(txy-tyx)*(txz-tzx)-(tyz-tzy)*(txx+tyy+tzz)) + 4*B/9*(txy*(5*txz-2*tzx)+5*tyx*(tzx+2*txz)+(5*tyz+tzy)*(tyy+tzz-2*txx))
J11 = 4*A/9*(-(txy-tyx)^2-(txz-tzx)^2+(tyz-tzy)^2+(txx+tyy+tzz)^2)+4*B/9*((txy-tyx)^2+(txz-tzx)^2-2(2tyz+tzy)*(tyz+2tzy)+2(txx+tyy-2tzz)*(txx-2tyy+tzz))
J33 = 4*A/9*(-(txz-tzx)^2-(txy-tyx)^2+(tyz-tzy)^2+(txx+tyy+tzz)^2)+4*B/9*((txz-tzx)^2+(tyz-tzy)^2-2(2txy+tyx)*(txy+2tyx)+2(tzz+tyy-2txx)*(txx-2tyy+tzz))
println("Jyz=",J23," Jxx=",J11," Jzz=",J33)
using Test
@test jex2[Symbol("σᵢʸ",",","σⱼᶻ")] ≈ J23
@test jex2[Symbol("σᵢˣ",",","σⱼˣ")] ≈ J11
@test jex2[Symbol("σᵢᶻ",",","σⱼᶻ")] ≈ J33


text:
const CS_ERROR_FORMAT_THRESHOLD = 0x00010000
const CS_FALSE = 0
const CS_SUCCESS = 1
const CS_CONFIG_CHANGED = 2
const CS_ASYNC_SUCCESS = 3
const CS_SEGMENTINFO_EMPTY = 4
const CS_NOT_INITIALIZED = -1
const CS_UNABLE_CREATE_RM = -3
const CS_INTERFACE_NOT_FOUND = -4
const CS_HANDLE_IN_USE = -5
const CS_INVALID_HANDLE = -6
const CS_INVALID_REQUEST = -7
const CS_NO_SYSTEMS_FOUND = -8
const CS_MEMORY_ERROR = -9
const CS_LOCK_SYSTEM_FAILED = -10
const CS_INVALID_STRUCT_SIZE = -11
const CS_INVALID_STATE = -12
const CS_INVALID_EVENT = -13
const CS_INVALID_SHARED_REGION = -14
const CS_INVALID_FILENAME = -15
const CS_SHARED_MAP_UNAVAILABLE = -16
const CS_INVALID_START = -17
const CS_INVALID_LENGTH = -18
const CS_SOCKET_NOT_FOUND = -19
const CS_SOCKET_ERROR = -20
const CS_NO_AVAILABLE_SYSTEM = -21
const CS_NULL_POINTER = -22
const CS_INVALID_CHANNEL = -23
const CS_INVALID_TRIGGER = -24
const CS_INVALID_EVENT_TYPE = -25
const CS_BUFFER_TOO_SMALL = -26
const CS_INVALID_PARAMETER = -27
const CS_INVALID_SAMPLE_RATE = -28
const CS_NO_EXT_CLK = -29
const CS_SEG_COUNT_TOO_BIG = -30
const CS_INVALID_SEGMENT_SIZE = -31
const CS_DEPTH_SIZE_TOO_BIG = -32
const CS_INVALID_CAL_MODE = -33
const CS_INVALID_TRIG_COND = -34
const CS_INVALID_TRIG_LEVEL = -35
const CS_INVALID_TRIG_SOURCE = -36
const CS_INVALID_EXT_TRIG = -37
const CS_INVALID_ACQ_MODE = -38
const CS_INVALID_IMPEDANCE = -39
const CS_INVALID_GAIN = -40
const CS_INVALID_COUPLING = -41
const CS_BUFFER_NOT_ALIGNED = -42
const CS_PRETRIG_DEPTH_TOO_BIG = -43
const CS_INVALID_TRIG_DEPTH = -44
const CS_FUNCTION_NOT_SUPPORTED = -45
const CS_HARDWARE_TIMEOUT = -46
const CS_INVALID_PARAMS_ID = -47
const CS_INVALID_POINTER_BUFFER = -48
const CS_CANNOT_LOCKDOWN_BUFFER = -49
const CS_DRIVER_ASYNC_REQUEST_BUSY = -50
const CS_INVALID_CHANNEL_COUNT = -51
const CS_INVALID_TRIGGER_COUNT = -52
const CS_INVALID_SEGMENT = -53
const CS_INVALID_SEGMENT_COUNT = -54
const CS_INVALID_CAPS_ID = -55
const CS_HANDLE_NOT_IN_USE = -56
const CS_INSUFFICIENT_RESOURCES = -57
const CS_INVALID_TRANSFER_MODE = -58
const CS_DRIVER_ASYNC_NOT_SUPPORTED = -59
const CS_INVALID_TRIGGER_ENABLED = -60
const CS_NOT_TRIGGER_FROM_SAME_CARD = -61
const CS_INVALID_PRETRIGGER_DEPTH = -62
const CS_INVALID_FW_VERSION = -63
const CS_INVALID_TRIGHOLDOFF = -64
const CS_INVALID_TRIGDELAY = -65
const CS_INVALID_STREAMING_BUFFER = -66
const CS_HW_FIFO_OVERFLOW = -67
const CS_INVALID_CARD = -68
const CS_INVALID_TOKEN = -69
const CS_MASTERSLAVE_DISCREPANCY = -70
const CS_INVALID_FIR_MULTIPLIER = -71
const CS_MINMAXDETECTQUEUE_INVALID = -72
const CS_ALLTRIGGERENGINES_USED = -73
const CS_MULREC_RAWDATA_TOOBIG = -74
const CS_INVALID_FILTER = -75
const CS_INVALID_POSITION = -76
const CS_EXT_CLK_OUT_OF_RANGE = -77
const CS_EXT_CLK_NOT_PRESENT = -78
const CS_RE_INIT_FAILED = -79
const CS_CHANNEL_PROTECT_FAULT = -80
const CS_INVALID_NUM_OF_AVERAGE = -81
const CS_INVALID_TIMESTAMP_CLOCK = -82
const CS_INVALID_FRM_CMD = -100
const CS_FRM_NO_RESPONSE = -101
const CS_INVALID_DAC_ADDR = -102
const CS_INVALID_EDGE = -103
const CS_INVALID_SELF_TEST = -104
const CS_NIOS_FAILED = -105
const CS_NO_INTERRUPT = -106
const CS_ADDONINIT_ERROR = -107
const CS_INVALID_ADC_ADDR = -108
const CS_ADC_ACCESS_ERROR = -109
const CS_DAC_CALIB_FAILURE = -200
const CS_CAL_BUSY_TIMEOUT = -201
const CS_TIMING_CAL_FAILED = -202
const CS_COARSE_OFFSET_CAL_FAILED = -203
const CS_FINE_OFFSET_CAL_FAILED = -204
const CS_GAIN_CAL_FAILED = -205
const CS_POSITION_CAL_FAILED = -206
const CS_CALIB_DAC_OUT_OF_RANGE = -207
const CS_CALIB_ADC_CAPTURE_FAILURE = -208
const CS_CALIB_MEAS_ADC_FAILURE = -209
const CS_CALIB_ADC_READ_FAILURE = -210
const CS_CHANNELS_NOTCALIBRATED = -211
const CS_MASTERSLAVE_CALIB_FAILURE = -212
const CS_EXTTRIG_CALIB_FAILURE = -213
const CS_INVALID_DACCALIBTABLE = -214
const CS_ADCALIGN_CALIB_FAILURE = -215
const CS_CALIB_REF_FAILURE = -216
const CS_FAST_CALIB_FAILURE = -217
const CS_ADC_DPA_FAILURE = -218
const CS_EXT_TRIG_DPA_FAILURE = -219
const CS_CLK_NOT_LOCKED = -220
const CS_INVALID_CALIB_BUFFER = -221
const CS_NULL_OFFSET_CAL_FAILED = -222
const CS_ADC_SKEW_CAL_FAILED = -223
const CS_ADC_PHASE_CAL_FAILED = -224
const CS_ACCAL_NOT_LOCKED = -225
const CS_NVRAM_NOT_INIT = -300
const CS_FLASH_NOT_INIT = -301
const CS_EEPROM_NOT_INIT = -302
const CS_ADDON_NOT_CONNECTED = -303
const CS_MS_BRIDGE_FAILED = -304
const CS_GIO_DETECT_FAILED = -305
const CS_ADDON_FPGA_LOAD_FAILED = -306
const CS_CSI_FILE_ERROR = -407
const CS_INVALID_DIRECTORY = -408
const CS_DISK_STREAM_NOT_INITIALIZED = -409
const CS_INVALID_DISK_STREAM_ACQ_COUNT = -410
const CS_INVALID_SIG_HEADER = -411
const CS_MISC_ERROR = -32767


text:
# Based on unbounded work-stealing deque:
# https://github.com/tkf/ConcurrentCollections.jl/blob/master/src/workstealing.jl

struct CircularVector{T} <: AbstractVector{T}
    log2length::Int
    data::Vector{T}
end

function CircularVector{T}(log2length::Int) where {T}
    @assert log2length >= 0
    data = Vector{T}(undef, 1 << log2length)
    T <: Signed && fill!(data, typemin(T))
    return CircularVector{T}(log2length, data)
end

Base.size(A::CircularVector{T}) where {T} = size(A.data)

@noinline function _check_data_length(A::CircularVector)
    @assert 1 << A.log2length == length(A.data)
end

Base.@propagate_inbounds function indexof(A::CircularVector, i::Int)
    @boundscheck _check_data_length(A)
    return (i - 1) & ((1 << A.log2length) - 1) + 1
end

Base.@propagate_inbounds Base.getindex(A::CircularVector, i::Int) = A.data[indexof(A, i)]

Base.@propagate_inbounds function Base.setindex!(A::CircularVector, v, i::Int)
    v = convert(eltype(A), v)
    A.data[indexof(A, i)] = v
end

mutable struct BoundedWorkStealingDeque{T}
    buffer::CircularVector{T}
    @atomic top::Int
    @atomic bottom::Int
    # TODO: pad
end

BoundedWorkStealingDeque{T}(size::Integer) where {T} =
    BoundedWorkStealingDeque{T}(CircularVector{T}(ceil(Int, log2(size))), 1, 1)

Base.eltype(::Type{BoundedWorkStealingDeque{T}}) where {T} = T

function Base.length(deque::BoundedWorkStealingDeque)
    bottom = @atomic deque.bottom
    top = @atomic deque.top
    return bottom - top
end

function isfull(deque::BoundedWorkStealingDeque)
    buffer = deque.buffer
    return length(deque) >= length(buffer) - 1
end

function trypush!(deque::BoundedWorkStealingDeque, v)
    v = convert(eltype(deque), v)
    bottom = @atomic deque.bottom
    top = @atomic deque.top
    buffer = deque.buffer
    current_size = bottom - top
    if current_size >= length(buffer) - 1
        return false
    end
    buffer[bottom] = v
    bottom += 1
    @atomic deque.bottom = bottom
    return true
end

function maybepop!(deque::BoundedWorkStealingDeque)
    bottom = @atomic deque.bottom
    buffer = deque.buffer
    bottom -= 1
    @atomic deque.bottom = bottom
    top = @atomic deque.top
    next_size = bottom - top
    if next_size < 0
        @atomic deque.bottom = top
        return nothing
    end
    r = Some(buffer[bottom])
    if next_size > 0
        return r
    end
    bottom = top + 1
    if !@atomicreplace(deque.top, top => top + 1)[2]
        r = nothing
    end
    @atomic deque.bottom = bottom
    return r
end

function maybepopfirst!(deque::BoundedWorkStealingDeque)
    top = @atomic deque.top
    bottom = @atomic deque.bottom
    buffer = deque.buffer
    current_size = bottom - top
    if current_size <= 0
        return nothing
    end
    r = Some(buffer[top])
    if @atomicreplace(deque.top, top => top + 1)[2]
        return r
    else
        return nothing
    end
end

function Base.pop!(deque::BoundedWorkStealingDeque)
    r = maybepop!(deque)
    if r === nothing
        error("deque is empty")
    else
        return something(r)
    end
end

function Base.popfirst!(deque::BoundedWorkStealingDeque)
    r = maybepopfirst!(deque)
    if r === nothing
        error("deque is empty")
    else
        return something(r)
    end
end


text:
macro timer(ex)
    quote
        local stats = Base.gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        local diff = Base.GC_Diff(Base.gc_num(), stats)
        timer_print(elapsedtime, diff.allocd, diff.total_time,
                   Base.gc_alloc_count(diff))
        val
    end
end

function timer_print(elapsedtime, bytes, gctime, allocs)
    @printf(STDERR, "%10.6f seconds", elapsedtime/1e9)
    if bytes != 0 || allocs != 0
        bytes, mb = Base.prettyprint_getunits(bytes, length(Base._mem_units), Int64(1024))
        allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000))
        if ma == 1
            @printf(STDERR, " (%d%s allocation%s: ", allocs, Base._cnt_units[ma], allocs==1 ? "" : "s")
        else
            @printf(STDERR, " (%.2f%s allocations: ", allocs, Base._cnt_units[ma])
        end
        if mb == 1
            @printf(STDERR, "%d %s%s", bytes, Base._mem_units[mb], bytes==1 ? "" : "s")
        else
            @printf(STDERR, "%.3f %s", bytes, Base._mem_units[mb])
        end
        if gctime > 0
            @printf(STDERR, ", %.2f%% gc time", 100*gctime/elapsedtime)
        end
        print(STDERR, ")")
    elseif gctime > 0
        @printf(STDERR, ", %.2f%% gc time", 100*gctime/elapsedtime)
    end
    println(STDERR)
end


text:
include("../utils.jl")

@testset "Essential manifold" begin
    M = EssentialManifold()
    a = π / 6
    b = π / 6
    c = π / 6
    r1 = [1.0 0.0 0.0; 0.0 cos(a) -sin(a); 0.0 sin(a) cos(a)]
    r2 = [cos(b) 0.0 sin(b); 0.0 1.0 0.0; -sin(b) 0.0 cos(b)]
    r3 = [cos(c) -sin(c) 0.0; sin(c) cos(c) 0.0; 0.0 0.0 1.0]
    nr = [1.0 0.0 0.0; 0.0 0.0 -1.0; 0.0 -1.0 0.0]
    p1 = [r1, r2]
    p2 = [r1, r3]
    p3 = [r2, r2]
    @testset "Essential manifold Basics" begin
        @test M.manifold == Rotations(3)
        @test repr(M) == "EssentialManifold(true)"
        @test manifold_dimension(M) == 5
        np1 = [r1, nr]
        np2 = [nr, nr]
        np3 = [r1, r2, r3]
        @test !is_point(M, r1)
        # first two components of r1 are not rotations
        @test_throws DomainError is_point(M, r1, true)
        @test_throws DomainError is_point(M, np3, true)
        @test is_point(M, p1)
        @test_throws ComponentManifoldError is_point(M, np1, true)
        @test_throws CompositeManifoldError is_point(M, np2, true)
        @test !is_vector(M, p1, 0.0)
        @test_throws DomainError is_vector(
            M,
            p1,
            [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],
            true,
        )
        @test !is_vector(M, np1, [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0])
        @test !is_vector(M, p1, p2)
        # projection test
        @test is_vector(M, p1, project(M, p1, log(M, p1, p2)))
        @test is_vector(M, p1, project(M, p1, log(M, p2, p1)))
    end
    @testset "Signed Essential" begin
        test_manifold(
            M,
            [p1, p2, p3],
            test_vector_spaces=true,
            test_project_point=true,
            projection_atol_multiplier=10,
            test_project_tangent=false, # since it includes vert_proj.
            test_musical_isomorphisms=false,
            test_default_vector_transport=true,
            test_representation_size=false,
            test_exp_log=true,
            mid_point12=nothing,
            exp_log_atol_multiplier=4,
            test_inplace=true,
            parallel_transport=true,
        )
    end
    @testset "Unsigned Essential" begin
        test_manifold(
            EssentialManifold(false),
            [p1, p2, p3],
            test_vector_spaces=true,
            test_project_point=true,
            projection_atol_multiplier=10,
            test_project_tangent=false, # since it includes vert_proj.
            test_musical_isomorphisms=false,
            test_default_vector_transport=true,
            test_representation_size=false,
            test_exp_log=true,
            mid_point12=nothing,
            exp_log_atol_multiplier=4,
            parallel_transport=true,
        )
    end
end


text:
using ITensors
using OptimKit
using Zygote
using Random: seed!

include("circuit.jl")

function gate(::Val{:Ry}; θ)
  return [
    cos(θ / 2) -sin(θ / 2)
    sin(θ / 2) cos(θ / 2)
  ]
end

function gate(::Val{:CX})
  return [
    1 0 0 0
    0 1 0 0
    0 0 0 1
    0 0 1 0
  ]
end

function Rylayer(N, θ⃗)
  return [("Ry", (n,), (θ=θ⃗[n],)) for n in 1:N]
end

function CXlayer(N)
  return [("CX", (n, n + 1)) for n in 1:2:(N - 1)]
end

# The variational circuit we want to optimize
function variational_circuit(θ⃗)
  N = length(θ⃗)
  return vcat(Rylayer(N, θ⃗), CXlayer(N), Rylayer(N, θ⃗), CXlayer(N))
end

N = 4
θ⃗ = 2π .* rand(N)
gates = variational_circuit(θ⃗)

s = siteinds("Qubit", N)
ψₘₚₛ = MPS(s, "0")
ψ = prod(ψₘₚₛ)
U = buildcircuit(gates, s)
# Create the target state
Uψ = apply(U, ψ)

@show inner_circuit(Uψ, U, ψ);

function loss(θ⃗)
  gates = variational_circuit(θ⃗)
  U = buildcircuit(gates, s)
  return -abs(inner_circuit(Uψ, U, ψ))^2
end

seed!(1234)
θ⃗₀ = randn!(copy(θ⃗))
@show θ⃗
@show loss(θ⃗)
@show loss'(θ⃗)
@show θ⃗₀
@show loss(θ⃗₀)
@show loss'(θ⃗₀)

fg(x) = (loss(x), convert(Vector, loss'(x)))
θ⃗ₒₚₜ, fₒₚₜ, gₒₚₜ, numfg, normgradhistory = optimize(fg, θ⃗₀, GradientDescent())
@show θ⃗ₒₚₜ
@show fₒₚₜ
@show gₒₚₜ
@show numfg
println("normgradhistory = ")
display(normgradhistory)

nothing


text:
module ImplicitBias

# Write your package code here.
export twonet

include("exp_bothlayers.jl")

end

max_stars_repo_name: brechetp/ImplicitBias

max_issues_repo_name: brechetp/ImplicitBias

max_forks_repo_name: brechetp/ImplicitBias


text:
"""
    indent(fun, io, n[; indent_first=true, kwargs...])

Print all the output of `fun` to `io` indented by `n` spaces. If
`!indent_first`, the first line is _not_ indented. `kwargs...` are
passed on to `printstyled`.

# Examples

```jldoctest
julia> indent(stdout, 6) do io
           println(io, "Hello")
           println(io, "World")
       end
      Hello
      World
```
"""
function indent(fun::Function, io::IO, n::Int; indent_first=true, kwargs...)
    indentation = repeat(" ", n)
    for (i,l) in enumerate(split(rstrip(redirect_output(fun)), "\n"))
        (indent_first || i > 1) && write(io, indentation)
        printstyled(io, l; kwargs...)
        println(io)
    end
end

"""
    indent(fun, io, first_line::String[; kwargs...])

Variant of [`indent`](@ref) that will first print `first_line` and
then print all the output from `fun` indented by the length of
`first_line`. `kwargs` are passed on to `printstyled`.

# Examples

```jldoctest
julia> indent(stdout, "Important information: ") do io
           println(io, "Hello")
           println(io, "World")
       end
Important information: Hello
                       World

julia> indent(stdout, "Important information: ") do io
           print_boxed(io) do io
               println(io, "Hello")
               println(io, "World")
           end
       end

Important information: ┌  Hello
                       └  World
```
"""
function indent(fun::Function, io::IO, first_line::String; kwargs...)
    printstyled(io, first_line; kwargs...)
    n = textwidth(last(split(first_line, "\n")))
    indent(fun, io, n; indent_first=false, kwargs...)
end


text:
using Plots; pyplot(fmt = :png)

include("ReadGeant4Hits.jl")
include("SetupSSD.jl")
include("EventSimulation.jl")
include("WaveformAnalysis.jl")
include("DigitalFilters.jl")
include("Plotting.jl")
include("SummaryReport.jl")

@info "Start"

temp = "140K"

# Setup the simulation
simFilename = "save/1Hex110K.h5"
geometryFile = "config_files/Si 1Hexagon 10.3-0.1.json"
chargeDriftConfigFile = string("config_files/drift_velocity_si_", temp, ".json")
#charge_drift_model = ADLChargeDriftModel(chargeDriftConfigFile)
CCDName = "config_files/1Hexagon.jl"
max_ref = 3

# Get either a saved simulation or a new one
sim = ReadSimulation(simFilename)
if isnothing(sim)
    sim = SetupDefaultSimulation(geometryFile, CCDName)
    CalculateDefaultDetectorFields!(sim, chargeDriftConfigFile, max_ref)
    SaveSimulation(simFilename, sim)
end

SetChargeDriftModel!(sim, chargeDriftConfigFile)
#SetChargeDriftModel!(sim, ADLChargeDriftModel())

geantDirectory = "/home/Group/NabSim/data/"
particle = "e-_40keV"
angles = [0,10,20,30,40,50,60,70,80]
geantFilenames = [string(geantDirectory, particle, "_", angle,"inc.root") for angle in angles]
gdfs = [GetHitInformation(file) for file in geantFilenames]

offset0 = CartesianPoint{T}(0, 0, 0)
events0 = []
for gdf in gdfs
    global events0 =[events0; DriftGeant4Events(gdf, sim, offset0, time_step=0.1u"ns",stepLimiter=1000000000)]
end
riseTimes0 = CollectRiseTimesFiltered(events0, dt=0.1e-9)

histogram(riseTimes0 / 1u"ns", bins = 50, alpha = 0.5, label=string(particle," ",temp))
savefig(string("plots/", particle, string("RiseTimeFiltered",temp,".png")))

@info "Stop"

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP