text:
#
# ```@setup setup
# # all examples are run on an example global.toml and downsized maps.
# ARGS = ["example.toml"]
# ```
#
# # [Setup (setup.jl)](@id setup)
#
# This pipeline is written in Julia, so you will need a [Julia](https://julialang.org/)
# installation in order to run the components. We recommend
# you use the precompiled binaries provided on the Julia website. Make sure to add the
# Julia executable to your path, as described in the
# [platform-specific instructions.](https://julialang.org/downloads/platform/)
#
#

# ## Package Installation
# We use the package manager in the Julia interpeter
# to install the latest versions of Healpix and PowerSpectra. This will be simpler in
# the future, when we tag a stable version of these packages for the General Registry.
# For now, we add the latest versions of these packages from GitHub. Note that package
# installation requires an internet connection, so unlike the other parts of the pipeline,
# `setup.jl` requires an internet connection. If you're on a cluster, that means you need
# to run this file on the head node in order to install packages.
using Pkg
Pkg.add.(["Healpix", "PowerSpectra", "CSV", "DataFrames", "TOML",
  "BlackBoxOptim", "FileIO", "JLD2", "FITSIO",
  "DataInterpolations", "Optim", "GaussianProcesses"]);
#
#
# The command-line interface for this basic pipeline setup script is
# ```@raw html
# <pre class="shell">
# <code class="language-shell hljs">$ julia setup.jl global.toml</code></pre>
# ```
# * It displays the contents of the global TOML configuration file named *global.toml*.
# * This script downloads the Planck data to the specified directories in *global.toml*.
#
#
#
# ## Configuration
# All of the pipeline scripts take a configuration TOML file as the first argument.
# We now print out just the `[dir]` entry in the TOML, which is what you will need to
# configure.
#

#src   This has the @example manually specified so that it runs. By default, nothing else
#src   in this script should run. The #src after makes it run in the script too.

# ```@example setup
# using TOML
# configfile = ARGS[1]  # read in the first command line argument
# println("config filename: ", configfile, "\n")
#
# # take a look at the config
# config = TOML.parsefile(configfile)
# TOML.print(Dict("dir"=>config))  # print just the "dir" TOML entry
# ```

##src
using TOML                                                                  #src
configfile = ARGS[1]  # read in the first command line argument             #src
println("config filename: ", configfile, "\n")                              #src
config = TOML.parsefile(configfile)                                         #src
TOML.print(config)  # print just the "dir" TOML entry   #src

# * The `scratch` space is where intermediary files are deposited.
# * Note that each map file has an identifier. This shortens the long names, but more
#   importantly allows one to set up a custom system of names when we cross-correlate Planck
#   with other experiments.
# * In this case, we preface all Planck maps and masks with `P`, and include the frequency
#   and split.

# # Downloading the Planck 2018 Data

## set up download urls
using Downloads

if length(ARGS) == 0
    TARGET_DIR = pwd()
else
    TARGET_DIR = ARGS[1]
end

## set up directories
mapdest = joinpath(config["scratch"], "maps")
maskdest = joinpath(config["scratch"], "masks")
beamdest = joinpath(config["scratch"], "beams")
mkpath(mapdest)
mkpath(maskdest)
mkpath(beamdest)

function download_if_necessary(url, dest; verbose=true)
    if isfile(dest) == false
        verbose && println("Downloading ", dest)
        @time Downloads.download(url, dest)
    else
        verbose && println("Extant, skip ", dest)
    end
end

#
## now read from config and then actually download
mapfiles = values(config["map"])
maskfiles = [values(config["maskT"])..., values(config["maskP"])...]

for f in mapfiles
    download_if_necessary(joinpath(config["url"]["maps"], f), joinpath(mapdest, f))
end

for f in maskfiles
    download_if_necessary(joinpath(config["url"]["masks"], f), joinpath(maskdest, f))
end

## just download beamfile to the target directory base.
beamfile = "HFI_RIMO_BEAMS_R3.01.tar.gz"
fullbeamfile = joinpath(beamdest, beamfile)
download_if_necessary(joinpath(config["url"]["beams"], beamfile), fullbeamfile)
run(`cp $(fullbeamfile) $(joinpath(beamdest, "tempbeamgz"))`)
run(`gzip -f -d $(joinpath(beamdest, beamfile))`)
run(`mv $(joinpath(beamdest, "tempbeamgz")) $(fullbeamfile)`)
beamfiletar = replace(beamfile, ".tar.gz"=>".tar")
run(`tar -xf $(joinpath(beamdest, beamfiletar)) --overwrite -C $(beamdest)`);

#
## We also want to retrieve the plic templates.
plicrefdir = joinpath(config["scratch"])
plicreffile = joinpath(config["scratch"], "plicref.tar.gz")
download_if_necessary(
    "https://github.com/xzackli/PSPipePlanckRender.jl/releases/download/0.1.2/plicref.tar.gz",
    plicreffile)
run(`tar -xzvf $(plicreffile) --overwrite -C $(plicrefdir)`);

# For the sims, we will use the full pixel weights to efficiently compute the spherical harmonic
# transforms. The pixel weights are lazily downloaded in the Healpix package. We trigger
# that downloading here. Similarly we trigger some downloading of example data from PowerSpectra.jl
# for the purpose of plotting.

using Healpix
Healpix.applyFullWeights!(HealpixMap{Float64, RingOrder}(2048))
using PowerSpectra
PowerSpectra.planck256_beamdir()

max_forks_repo_name: sgiardie/PSpipe


text:
# ---
# title: 841. Keys and Rooms
# id: problem841
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Graph
# link: <https://leetcode.com/problems/keys-and-rooms/description/>
# hidden: true
# ---
#
# There are `N` rooms and you start in room `0`.  Each room has a distinct
# number in `0, 1, 2, ..., N-1`, and each room may have some keys to access the
# next room.
#
# Formally, each room `i` has a list of keys `rooms[i]`, and each key
# `rooms[i][j]` is an integer in `[0, 1, ..., N-1]` where `N = rooms.length`.  A
# key `rooms[i][j] = v` opens the room with number `v`.
#
# Initially, all the rooms start locked (except for room `0`).
#
# You can walk back and forth between rooms freely.
#
# Return `true` if and only if you can enter every room.
#
# **Example 1:**
#
#
#
#     Input: [[1],[2],[3],[]]
#     Output: true
#     Explanation:
#     We start in room 0, and pick up key 1.
#     We then go to room 1, and pick up key 2.
#     We then go to room 2, and pick up key 3.
#     We then go to room 3.  Since we were able to go to every room, we return true.
#
#
# **Example 2:**
#
#
#
#     Input: [[1,3],[3,0,1],[2],[0]]
#     Output: false
#     Explanation: We can't enter the room with number 2.
#
#
# **Note:**
#
#   1. `1 <= rooms.length <= 1000`
#   2. `0 <= rooms[i].length <= 1000`
#   3. The number of keys in all rooms combined is at most `3000`.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
import Weave
using DataFrames
using JSON
using RCall
import PlotlyJS

function jhtml(source::String)
    Weave.weave(source,doctype="pandoc",cache=:user)
    ofile = split(basename(source),".")[1] * ".md"
    htmlfile = split(basename(source),".")[1] * ".html"
    pdir = dirname(@__FILE__)
    vtemplate = joinpath(pdir, "../templates/vuetify.html")
    vcss = read(joinpath(pdir, "../templates/vuetify.css"),String)
    prismjs = read(joinpath(pdir, "../templates/prism.js"),String)
    pandoc_result = run(`pandoc -f markdown -t html $ofile -o $htmlfile --katex --section-divs --template $vtemplate --toc --css $vcss --no-highlight -V prismjs=$prismjs` );
    return nothing;
end

function htmltable(df;title="",cnames = names(df))
  function df2json(df::DataFrame)
    [Dict([string(index) => (df[!, index][i]) for index in names(df)]) for i = 1:nrow(df)] |> JSON.json
  end
  trows = df2json(df)
  thead = JSON.json([Dict("value"=>names(df)[i],"text"=>cnames[i]) for i in 1:length(names(df))] )
  temp_string = ""
  for i in names(df)
    temp_string *= """<template v-slot:item.$i="{ item }"></template>
    """
  end
  return """
  ```{=html}
  <vue-data-table>
  <template slot="title"><span style="color:#4e6b95; font-weight:normal">$title</span></template>
  <template slot-scope="{search, cfilter}">
  <v-data-table dense
      :headers=\'$thead\'
      :items=\'$trows\'
      :search = "search"
      :custom-filter = "cfilter"
      multi-sort
      class="elevation-1"
    >
    $temp_string
    </v-data-table>
  </template>
  </vue-data-table>
  ```"""
end
function jggplot(rplot,df,fname::String;type::String = "markdown", width = 8, height =5)
    jggplot(rplot,df; width = width, height = height)
    R"""
    if (!dir.exists("figures")) {dir.create("figures")}
    ggsave(filename = paste0("figures/",$fname), width = $width, height = $height,type = "cairo", dpi = 300)
    """
    if type == "markdown"
      return "![](figures/$fname)"
    else
      return ""
    end
end

function jggplot(rplot,df;width = 8, height =5)
  @rput(df)
  R"library(hrbrthemes)"
  R"library(ggplot2)"
  p = reval(rplot)
end

function plotlygraph(a,pid)
  pdata = json(a.plot)
  println("""```{=html}
    <plotly-vue pid=$pid pdata = '$pdata'></plotly-vue>
  ```""")
end


text:
function ccf_orders(order_list_timeseries::AbstractChunkListTimeseries,  line_list_df::DataFrame, pipeline::PipelinePlan; verbose::Bool = false, calc_ccf_var::Bool = false, recalc::Bool = false,
   orders_to_use = RvSpectMLBase.orders_to_use_default(order_list_timeseries.inst),
   range_no_mask_change::Real=RvSpectMLBase.max_bc, ccf_mid_velocity::Real=0.0,
    v_step::Real=250.0, v_max=RvSpectMLBase.max_bc, mask_scale_factor::Real=1,
    mask_type::Symbol = :tophat, Δfwhm::AbstractVector{T} = zeros(0),
    allow_nans::Bool = false ) where { T<:Real }

   if need_to(pipeline, :ccf_orders)  || recalc # Compute order CCF's & measure RVs
      if mask_type == :tophat
         mask_shape = CCFs.TopHatCCFMask(order_list_timeseries.inst, scale_factor=mask_scale_factor)
      elseif mask_type == :gaussian
         mask_shape = CCFs.GaussianCCFMask(order_list_timeseries.inst, σ_scale_factor=mask_scale_factor)
      elseif mask_type == :supergaussian
         mask_shape = CCFs.SuperGaussianCCFMask(order_list_timeseries.inst, σ_scale_factor=mask_scale_factor)
      elseif mask_type == :halfcos
         mask_shape = CCFs.CosCCFMask(order_list_timeseries.inst, scale_factor=mask_scale_factor)
      else
         @error("Requested mask shape (" * string(mask_type) * " not avaliable.")
      end

      line_list = hasproperty(line_list_df, :order) ? CCFs.BasicLineList2D(line_list_df.lambda, line_list_df.weight, line_list_df.order) :
                  CCFs.BasicLineList(line_list_df.lambda, line_list_df.weight )
      ccf_plan = BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=range_no_mask_change, step=v_step, max=v_max, allow_nans=allow_nans)
      v_grid = calc_ccf_v_grid(ccf_plan)
      tstart = now()    # Compute CCFs for each order

      if calc_ccf_var
         (order_ccfs, order_ccf_vars ) = EchelleCCFs.calc_order_ccf_and_var_chunklist_timeseries(order_list_timeseries, ccf_plan, Δfwhm=Δfwhm)
         set_cache!(pipeline,:ccf_orders, (order_ccfs = order_ccfs, order_ccf_vars = order_ccf_vars, v_grid=v_grid) )
      else
         order_ccfs = EchelleCCFs.calc_order_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan, Δfwhm=Δfwhm)
         set_cache!(pipeline,:ccf_orders, (order_ccfs = order_ccfs, v_grid=v_grid) )
      end
      if verbose   println("# Order CCFs runtime: ", now()-tstart)  end

      if save_data(pipeline, :ccf_orders)
         for (i, order) in orders_to_use
            if !(sum(order_ccfs[:,i,:]) > 0)   continue    end
            local t = Tables.table( order_ccfs[:,i,:]', header=Symbol.(v_grid) )
            CSV.write(joinpath(output_dir,target_subdir * "_ccf_order=" * string(order) * ".csv"),t)
            if calc_ccf_var
               local t = Tables.table( order_ccf_vars[:,i,:]', header=Symbol.(v_grid) )
               CSV.write(joinpath(output_dir,target_subdir * "_ccf_var_order=" * string(order) * ".csv"),t)
            end
         end
      end
      dont_need_to!(pipeline, :ccf_orders);
   end
   if has_cache(pipeline,:ccf_orders) return read_cache(pipeline,:ccf_orders)
   else   @error("Invalid pipeline state.")          end
end


text:
using GeometricIntegrators.Config
using GeometricIntegrators.Solvers
using Test

include("linear_solvers_tests.jl")
include("jacobian_tests.jl")
include("nonlinear_solvers_tests.jl")


text:
AbstractJordan = Union{AbstractClosedCurve,AbstractClosedPath}
abstract type AbstractRegion end

# Required methods
boundary(R::AbstractRegion) = @error "No boundary() method defined for type $(typeof(R))"

"""
	in(z::Number,R::AbstractRegion;tol=<default>)
	z ∈ R   (type "\\in" followed by tab)
True if `z` is in the region `R`.
"""
in(z::Number,R::AbstractRegion;tol=DEFAULT[:tol]) = @error "No in() method defined for type $(typeof(R))"

"""
	isfinite(R::AbstractRegion)
Return `true` if the region is bounded in the complex plane.
"""
isfinite(R::AbstractRegion) = @error "No isfinite() method defined for type $(typeof(R))"

# Default implementations

"""
	(type) RegionIntersection
Representation of the intersection of two regions.
"""
struct RegionIntersection <: AbstractRegion
	one::AbstractRegion
	two::AbstractRegion
end
in(z::Number,R::RegionIntersection) = in(z,R.one) && in(z,R.two)

"""
	(type) RegionUnion
Representation of the union of two regions.
"""
struct RegionUnion <: AbstractRegion
	one::AbstractRegion
	two::AbstractRegion
end
in(z::Number,R::RegionUnion) = in(z,R.one) || in(z,R.two)

"""
	intersect(R1::AbstractRegion,R2::AbstractRegion)
	R1 ∩ R2    (type "\\cap" followed by tab key)
Create the region that is the intersection of `R1` and `R2`.
"""
intersect(R1::AbstractRegion,R2::AbstractRegion) = RegionIntersection(R1,R2)

"""
	union(R1::AbstractRegion,R2::AbstractRegion)
	R1 ∪ R2    (type "\\cup" followed by tab key)
Create the region that is the union of `R1` and `R2`.
"""
union(R1::AbstractRegion,R2::AbstractRegion) = RegionUnion(R1,R2)

#
# AbstractConnectedRegion
#

abstract type AbstractConnectedRegion{N} <: AbstractRegion end

# Required methods
innerboundary(R::AbstractConnectedRegion) = @error "No innerboundary() method defined for type $(typeof(R))"
outerboundary(R::AbstractConnectedRegion) = @error "No outerboundary() method defined for type $(typeof(R))"
boundary(R::AbstractConnectedRegion) = outerboundary(R),innerboundary(R)

# Default implementations

+(R::AbstractConnectedRegion,z::Number) = typeof(R)(outerboundary(R)+z,innerboundary(R).+z)
+(z::Number,R::AbstractConnectedRegion) = +(R,z)

-(R::AbstractConnectedRegion) = typeof(R)(-outerboundary(R),-innerboundary(R))
-(R::AbstractConnectedRegion,z::Number) = +(R,-z)
-(z::Number,R::AbstractConnectedRegion) = +(z,-R)

*(R::AbstractConnectedRegion,z::Number) = typeof(R)(outerboundary(R)*z,innerboundary(R)*z)
*(z::Number,R::AbstractConnectedRegion) = R*z

/(R::AbstractConnectedRegion,z::Number) = *(R,1/z)
#/(z::Number,R::AbstractConnectedRegion) = z*inv(R)
#inv(p::AbstractConnectedRegion) = typeof(R)([inv(c) for c in curves(p)])

#
# concrete implementations
#

#
# ExteriorRegion
#
struct ExteriorRegion{N} <: AbstractConnectedRegion{N}
	inner::AbstractVector
	function ExteriorRegion{N}(inner) where N
		@assert N == length(inner) "Incorrect connectivity"
		@assert all(c isa AbstractJordan for c in inner) "Boundary components must be closed curves or paths"
		@assert all(isfinite.(inner)) "Inner boundaries must be finite"
		# correct orientations of inner components
		b = copy(inner)
		for c in b
			if isoutside(Inf,c)
				 c = reverse(c)
			end
		end
		new(b)
	end
end

ExteriorRegion(inner) = ExteriorRegion{length(inner)}(inner)
in(z::Number,R::ExteriorRegion) = all( isoutside(z,c) for c in R.inner )
innerboundary(R::ExteriorRegion) = R.inner
outerboundary(R::ExteriorRegion) = nothing

#
# ConnectedRegion
#

"""
	(type) ConnectedRegion{N}
Representation of a `N`-connected region in the extended complex plane.
"""
struct ConnectedRegion{N} <: AbstractConnectedRegion{N}
	outer::Union{Nothing,AbstractJordan}
	inner::AbstractVector
	function ConnectedRegion{N}(outer,inner) where N
		n = length(inner) + !isnothing(outer)
		@assert N == n "Incorrect connectivity"
		if !isnothing(outer)
			# correct orientation of outer component?
			isin = [isinside(point(c,0),outer) for c in inner ]
			if all(.!isin)
				outer = reverse(outer)
			else
				@assert all(isin) "Inner components appear to be crossing the outer boundary"
			end
		end
		# correct orientations of inner components?
		@assert all(isfinite.(inner)) "Inner boundaries must be finite"
		for c in inner
			if !isfinite(interior(c))
				c = reverse(c)
			end
		end
		new(outer,inner)
	end
end

"""
	ConnectedRegion(outer,inner)
Construct an open connected region by specifying its boundary components. The `outer` boundary could be `nothing` or a closed curve or path. The `inner` boundary should be a vector of one or more nonintersecting closed curves or paths. The defined region is interior to the outer boundary and exterior to all the components of the inner boundary, regardless of the orientations of the given curves.
"""
function ConnectedRegion(outer,inner)
	n = length(inner)
	if isnothing(outer) || isempty(outer)
		ExteriorRegion{n}(inner)
	else
		ConnectedRegion{n+1}(outer,inner)
	end
end

function in(z::Number,R::ConnectedRegion)
	all( isoutside(z,c) for c in R.inner ) && isinside(z,R.outer)
end

outerboundary(R::ConnectedRegion) = R.outer
innerboundary(R::ConnectedRegion) = R.inner
#innerboundary(R::ConnectedRegion{2}) = R.inner[1]

#
# special cases
#

"""
	between(outer,inner)
Construct the region interior to the closed curve or path `outer` and interior to `inner`.
"""
function between(outer::AbstractJordan,inner::AbstractJordan)
	if isfinite(outer) && isinside(Inf,outer)
		outer = reverse(outer)
	end
	if isfinite(inner) && isoutside(Inf,inner)
		inner = reverse(inner)
	end
	ConnectedRegion{2}(outer,[inner])
end

# Annulus

"""
	(type) Annulus
Representation of the region between two circles.
"""
struct Annulus <: AbstractConnectedRegion{2}
	outer::Circle
	inner::Circle
	function Annulus(outer::Circle,inner::Circle)
		@assert(outer.center ≈ inner.center)
		if isinside(Inf,outer)
			outer = reverse(outer)
		end
		if isoutside(Inf,inner)
			inner = reverse(inner)
		end
		new(outer,inner)
	end
end
"""
	Annulus(radouter,radinner)
	Annulus(radouter,radinner,center)
Construct a concentric annulus of outer radius `radouter` and inner radius `radinner` centered at `center`. If the center is not given, the origin is used.
"""
function Annulus(outerrad::Real,innerrad::Real,center::Number=0)
	@assert outerrad > innerrad > 0
	Annulus(Circle(center,outerrad,true),Circle(center,innerrad,false))
end

innerboundary(A::Annulus) = A.inner
outerboundary(A::Annulus) = A.outer
in(z::Number,A::Annulus) = isinside(z,A.outer) && isoutside(z,A.inner)
isfinite(::Annulus) = true

function show(io::IO,::MIME"text/plain",R::Annulus)
	print(io,"Annulus in the complex plane:\n")
	print(io,"   centered at ",R.outer.center," with distances from ",R.inner.radius," to ",R.outer.radius)
end


text:
abstract type SelectionHeuristic end

## HELPERS

measure_adjusted_weights(weights, measures) =
    if weights isa Nothing
        vcat([signature(measures[1]), ], zeros(length(measures) - 1))
    else
        length(weights) == length(measures) ||
            throw(DimensionMismatch(
                "`OptimizeAggregatedMeasurement` heuristic "*
                "is being applied to a list of measures whose length "*
                "differs from that of the specified `weights`. "))
        signature.(measures) .* weights
    end

## OPTIMIZE AGGREGATED MEASURE

"""
    NaiveSelection(; weights=nothing)

Construct a common selection heuristic for use with `TunedModel` instances
which only considers measurements aggregated over all samples (folds)
in resampling.

For each entry in the tuning history, one defines a penalty equal to
the evaluations of the `measure` specified in the `TunedModel`
instance, aggregated over all samples, and multiplied by `-1` if `measure`
is a `:score`, and `+`` if it is a loss. The heuristic declares as
"best" (optimal) the model whose corresponding entry has the lowest
penalty.

If `measure` is a vector, then the first element is used, unless
per-measure `weights` are explicitly specified. Weights associated
with measures that are neither `:loss` nor `:score` are reset to zero.

"""
struct NaiveSelection <: SelectionHeuristic
    weights::Union{Nothing, Vector{Real}}
end
NaiveSelection(; weights=nothing) =
    NaiveSelection(weights)

function best(heuristic::NaiveSelection, history)
    first_entry = history[1]
    measures = first_entry.measure
    weights = measure_adjusted_weights(heuristic.weights, measures)
    measurements = [weights'*(h.measurement) for h in history]
    measure = first(history).measure[1]
    if orientation(measure) == :score
                               measurements = -measurements

    end
    best_index = argmin(measurements)
    return history[best_index]
end

MLJTuning.supports_heuristic(::Any, ::NaiveSelection) =
    true


text:
using DataDeps
using CSV, DataFrames

register(DataDep(
    "TCGA-PANCAN-HiSeq",
    """
    Dataset: TCGA-PANCAN-HiSeq
    Author: Samuele Fiorini, University of Genoa
    Website: https://archive.ics.uci.edu/ml/datasets/gene+expression+cancer+RNA-Seq

    Redistributed under Creative Commons license (http://creativecommons.org/licenses/by/3.0/legalcode)

    Observations: 801
    Features:     20531
    Classes:      5
    """,
    "https://archive.ics.uci.edu/ml/machine-learning-databases/00401/TCGA-PANCAN-HiSeq-801x20531.tar.gz",
    "e6ea3628f9656efa2d7cf382bf5a403e2cc214df3d90312b1928f5650b43a559",
    post_fetch_method = (path -> begin
        DataDeps.unpack(path, keep_originals=false)

        # Read the data and labels files.
        tmpdir = "TCGA-PANCAN-HiSeq-801x20531"
        data = CSV.read(joinpath(tmpdir, "data.csv"), DataFrame, header=true)
        labels = CSV.read(joinpath(tmpdir, "labels.csv"), DataFrame, header=true)

        # Create temporary version with columns Column1, Class, gene_0, gene_1, ...
        tmpdf = innerjoin(labels, data, on=:Column1)

        # Process the file as usual.
        MVDA.process_dataset(tmpdf,
            target_index=2,
            feature_indices=3:ncol(tmpdf),
            ext=".csv.gz"
        )

        # Clean up.
        rm(tmpdir, recursive=true)
    end),
))


text:
"""Plotting functions for lattices and spins on lattices.
"""

function plot_lattice!(ax, lattice::Lattice; colors=:Set1_9, markersize=200, linecolor=:grey, linewidth=1.0, kwargs...)
    unique_types = unique(lattice.types)
    colors = GLMakie.to_colormap(colors, 9)

    # Plot markers at each site
    pts = GLMakie.Point3f.(vec(lattice))
    for (i, type) in enumerate(unique_types)
        basis_idxs = findall(isequal(type), lattice.types)
        GLMakie.scatter!(ax, pts; label=type, color=colors[i], markersize=markersize, kwargs...)
    end

    # For some odd reason, the sites will not appear unless this happens afterwards
    # Plot the unit cell mesh
    plot_cells!(ax, lattice; color=linecolor, linewidth=linewidth)
end

function _setup_scene(; show_axis=false)
    fig = GLMakie.Figure()
    ax = GLMakie.LScene(fig[1, 1], show_axis=show_axis)
    return fig, ax
end

function plot_lattice(lattice::Lattice; kwargs...)
    fig, ax = _setup_scene()
    plot_lattice!(ax, lattice; kwargs...)
    # TODO: Markers are often way too big.
    fig[1, 2] = GLMakie.Legend(fig, ax, "Species")
    fig
end

"""
    plot_lattice(crystal, latsize=(3,3,3); kwargs...)

Plots a crystal lattice with `latsize` unit cells along each
lattice vector. Other keyword arguments are:

colors=:Set1_9, markersize=20, linecolor=:grey, linewidth=1.0, kwargs...
# Arguments
- `linecolor=:grey`: Sets the colors on the unit cell guide lines
- `linewidth=1.0`  : Sets the width of the unit cell guide lines
- `markersize=20`  : Sets the size of the atomic sites
- `colors=:Set1_9` : Sets the colors used for the atomic sites

Additional keyword arguments are given to `GLMakie.scatter!` which
draws the points.
"""
plot_lattice(cryst::Crystal, latsize=(3,3,3); kwargs...) = plot_lattice(Lattice(cryst, latsize); kwargs...)

function plot_bonds(lattice::Lattice, ints::Vector{<:AbstractInteractionCPU};
                    colors=:Dark2_8, bondwidth=4, kwargs...)
    fig, ax = _setup_scene()

    # Plot the bonds, all relative to a central atom on the first sublattice
    # TODO: Make selectable in GUI
    basis_idx = 1

    colors = GLMakie.to_colormap(colors, 8)
    # Sort interactions so that longer bonds are plotted first
    sorted_ints = sort(
        ints,
        by=int->(
            iszero(length(int.bondtable))
            ? Inf
            : distance(lattice, first(all_bonds(int.bondtable))[1])
        )
    )

    # Plot the lattice
    plot_lattice!(ax, lattice; kwargs...)

    toggles = Vector{GLMakie.Toggle}()
    labels = Vector{GLMakie.Label}()
    cent_cell = CartesianIndex(div.(lattice.size .+ 1, 2)...)
    cent_pt = lattice[basis_idx, cent_cell]
    for (n, int) in enumerate(sorted_ints)
        if !isa(int, AbstractPairIntCPU)
            continue
        end

        pts = Vector{GLMakie.Point3f}()
        for (bond, _) in sublat_bonds(int.bondtable, basis_idx)
            new_cell = offset(cent_cell, bond.n, lattice.size)
            bond_pt = lattice[bond.j, new_cell]
            push!(pts, GLMakie.Point3f(cent_pt))
            push!(pts, GLMakie.Point3f(bond_pt))
        end
        if length(pts) == 0
            continue
        end

        color = colors[mod1(n, 8)]
        seg = GLMakie.linesegments!(pts; linewidth=bondwidth, label=int.label, color=color)
        tog = GLMakie.Toggle(fig, active=true)
        GLMakie.connect!(seg.visible, tog.active)
        push!(toggles, tog)
        push!(labels, GLMakie.Label(fig, int.label))
    end
    GLMakie.axislegend()
    if length(toggles) > 0
        fig[1, 2] = GLMakie.grid!(hcat(toggles, labels), tellheight=false)
    end
    fig
end

"""
    plot_bonds(cryst::Crystal, ints, latsize=(3,3,3), sites_info=SiteInfo[]; kwargs...)

Plot a list of pair interactions defined on a `Crystal`. `latsize` sets how many
unit cells are plotted, and `kwargs` are passed to to `plot_lattice!`.
"""
function plot_bonds(cryst::Crystal, ints::Vector{<:AbstractInteraction}, latsize=(3,3,3);
                    kwargs...)
    latsize = collect(Int64.(latsize))
    lattice = Lattice(cryst, latsize)
    all_sites_info = propagate_site_info(cryst, SiteInfo[])
    ℋ = HamiltonianCPU(ints, cryst, latsize, all_sites_info)
    pair_ints = vcat(ℋ.heisenbergs, ℋ.diag_coups, ℋ.gen_coups)
    plot_bonds(lattice, pair_ints; kwargs...)
end

"""
    plot_bonds(sys::SpinSystem; kwargs...)

Plot all pair interactions appearing in `sys.hamiltonian`, on the
underlying crystal lattice. `kwargs` are passed to `plot_lattice!`.
"""
@inline function plot_bonds(sys::SpinSystem; kwargs...)
    ℋ = sys.hamiltonian
    pair_ints = vcat(ℋ.heisenbergs, ℋ.diag_coups, ℋ.gen_coups)
    plot_bonds(sys.lattice, pair_ints; kwargs...)
end

"""
    plot_all_bonds(crystal::Crystal, max_dist, latsize=(3,3,3); kwargs...)

Plot all bond equivalency classes present in `crystal` up to a maximum
bond length of `max_dist`. `latsize` controls how many unit cells are
plotted along each axis. `kwargs` are passed to `plot_bonds`.
"""
function plot_all_bonds(crystal::Crystal, max_dist, latsize=(3,3,3); kwargs...)
    ref_bonds = reference_bonds(crystal, max_dist)
    interactions = QuadraticInteraction[]

    prev_dist = 0.0
    dist = 0
    class = 1
    for bond in ref_bonds
        # Exclude self (on-site) "bonds"
        if !(bond.i == bond.j && all(isequal(0), bond.n))
            if distance(crystal, bond) ≈ prev_dist
                class += 1
            else
                dist += 1
                class = 1
            end
            label = "J$(dist)_$(class)"
            push!(interactions, heisenberg(1.0, bond, label))
        end
    end
    @assert length(interactions) > 0 "No non-self interactions found!"

    plot_bonds(crystal, interactions, latsize; kwargs...)
end

"Plot all bonds between equivalent sites i and j"
function plot_all_bonds_between(crystal, i, j, max_dist, latsize=(3,3,3); kwargs...)
    ref_bonds = reference_bonds(crystal, max_dist)
    interactions = QuadraticInteraction[]

    prev_dist = 0.0
    dist = 0
    class = 1
    for bond in ref_bonds
        # Exclude self (on-site) "bonds"
        onsite = bond.i == bond.j && iszero(bond.n)
        target = bond.i == i && bond.j == j
        if !onsite && target
            if distance(crystal, bond) ≈ prev_dist
                class += 1
            else
                dist += 1
                class = 1
            end
            label = "J$(dist)_$(class)"
            push!(interactions, heisenberg(1.0, bond, label))
        end
    end
    @assert length(interactions) > 0 "No non-self interactions found!"

    plot_bonds(crystal, interactions, latsize; kwargs...)
end

"Plot the outlines of the unit cells of a lattice"
function plot_cells!(ax, lattice::Lattice; color=:grey, linewidth=1.0, kwargs...)
    lattice = brav_lattice(lattice)

    pts = Vector{GLMakie.Point3f}()
    nx, ny, nz = lattice.size
    for j in 1:ny
        for k in 1:nz
            bot_pt, top_pt = lattice[1, 1, j, k], lattice[1, nx, j, k]
            push!(pts, GLMakie.Point3f(bot_pt))
            push!(pts, GLMakie.Point3f(top_pt))
        end
        for i in 1:nx
            left_pt, right_pt = lattice[1, i, j, 1], lattice[1, i, j, nz]
            push!(pts, GLMakie.Point3f(left_pt))
            push!(pts, GLMakie.Point3f(right_pt))
        end
    end
    for k in 1:nz
        for i in 1:nx
            left_pt, right_pt = lattice[1, i, 1, k], lattice[1, i, ny, k]
            push!(pts, GLMakie.Point3f(left_pt))
            push!(pts, GLMakie.Point3f(right_pt))
        end
    end

    GLMakie.linesegments!(ax, pts; color=color, linewidth=linewidth)
end

function plot_spins(lat::Lattice, spins; linecolor=:grey, arrowcolor=:red,
                    linewidth=0.1, arrowsize=0.2, arrowlength=0.2, kwargs...)
    fig, ax = _setup_scene()

    pts = GLMakie.Point3f.(vec(lat))
    vecs = GLMakie.Vec3f.(vec(spins))

    GLMakie.arrows!(
        ax, pts, vecs;
        linecolor=linecolor, arrowcolor=arrowcolor, linewidth=linewidth, arrowsize=arrowsize,
        lengthscale=arrowlength, kwargs...
    )
    fig
end

"""
    plot_spins(sys::SpinSystem; linecolor=:grey, arrowcolor=:red, linewidth=0.1,
                                arrowsize=0.3, arrowlength=1.0, kwargs...)

Plot the spin configuration defined by `sys`. `kwargs` are passed to `GLMakie.arrows`.
"""
plot_spins(sys::SpinSystem; kwargs...) = plot_spins(sys.lattice, sys.sites; kwargs...)

# No support for higher than 3D visualization, sorry!

"""
    anim_integration(sys, fname, steps_per_frame, Δt, nframes; kwargs...)

Produce an animation of constant-energy Landau-Lifshitz dynamics of the given `sys`.
# Arguments:
- `sys::SpinSystem`: The spin system to integrate.
- `fname::String`: The path to save the animation to.
- `steps_per_frame::Int`: The number of integration steps to take per frame.
- `Δt::Float64`: The integration timestep size.
- `nframes::Int`: The number of frames to produce in the animation.

Other keyword arguments are passed to `GLMakie.arrows`.
"""
function anim_integration(
    sys::SpinSystem, fname, steps_per_frame, Δt, nframes;
    linecolor=:grey, arrowcolor=:red, linewidth=0.1, arrowsize=0.2, arrowlength=0.2,
    kwargs...
)
    fig, ax = _setup_scene()

    pts = GLMakie.Point3f.(vec(sys.lattice))
    vecs = GLMakie.Observable(GLMakie.Vec3f.(vec(sys)))
    GLMakie.arrows!(
        ax, pts, vecs;
        linecolor=linecolor, arrowcolor=arrowcolor, linewidth=linewidth, arrowsize=arrowsize,
        lengthscale=arrowlength, kwargs...
    )
    display(fig)

    framerate = 30
    integrator = HeunP(sys)

    GLMakie.record(fig, fname, 1:nframes; framerate=framerate) do frame
        for step in 1:steps_per_frame
            evolve!(integrator, Δt)
        end
        vecs[] = GLMakie.Vec3f.(vec(sys))
    end
end

"""
    live_integration(sys, steps_per_frame, Δt; kwargs...)

Performs endless live constant-energy Landau-Lifshitz integration
in an interactive window.
"""
function live_integration(
    sys::SpinSystem, steps_per_frame, Δt;
    linecolor=:grey, arrowcolor=:red, linewidth=0.1, arrowsize=0.2,
    arrowlength=0.2, framerate=30, kwargs...
)
    fig, ax = _setup_scene()

    pts = GLMakie.Point3f.(vec(sys.lattice))
    vecs = GLMakie.Observable(GLMakie.Vec3f.(vec(sys)))
    GLMakie.arrows!(
        ax, pts, vecs;
        linecolor=linecolor, arrowcolor=arrowcolor, linewidth=linewidth, arrowsize=arrowsize,
        lengthscale=arrowlength, kwargs...
    )
    display(fig)

    integrator = HeunP(sys)

    while true
        for step in 1:steps_per_frame
            evolve!(integrator, Δt)
        end
        vecs[] = GLMakie.Vec3f.(vec(sys))
        sleep(1/framerate)
    end
end

"""
    live_langevin_integration(sys, steps_per_frame, Δt, kT; α=0.1, kwargs...)

Performs endless live Langevin Landau-Lifshitz integration
in an interactive window.
"""
function live_langevin_integration(
    sys::SpinSystem, steps_per_frame, Δt, kT;
    linecolor=:grey, arrowcolor=:red, linewidth=0.1, arrowsize=0.2,
    arrowlength=0.2, α=0.1, framerate=30, kwargs...
)
    fig, ax = _setup_scene()
    pts = GLMakie.Point3f.(vec(sys.lattice))
    vecs = GLMakie.Observable(GLMakie.Vec3f.(vec(sys)))

    GLMakie.arrows!(
        ax, pts, vecs;
        linecolor=linecolor, arrowcolor=arrowcolor, linewidth=linewidth, arrowsize=arrowsize,
        lengthscale=arrowlength, kwargs...
    )
    display(fig)

    integrator = LangevinHeunP(sys, kT, α)

    while true
        for step in 1:steps_per_frame
            evolve!(integrator, Δt)
        end
        vecs[] = GLMakie.Vec3f.(vec(sys))
        sleep(1/framerate)
    end
end

"Plots slices of a 3D structure factor. Input array should have shape [3, Lx, Ly, Lz, T]"
function plot_3d_structure_factor(sfactor::Array{Float64, 5}, iz)
    fig, ax = _setup_scene(; show_axis=true)

    Sdim, Lx, Ly, Lz, T = size(sfactor)
    # Average over Sxx, Syy, Szz - in future give controls to user
    sfactor = dropdims(sum(sfactor, dims=1) ./ 3, dims=1)
    # Index out the asked-for slice - in future give controls to user
    sfactor = sfactor[:, :, iz, :]

    # Nearest-neighbor upsample this to a similar size as the T axis.
    sampx, sampy = div(T, 2Lx), div(T, 2Ly)
    sfactor = repeat(sfactor, inner=(sampx, sampy, 1))

    kx = 1:sampx*Lx
    ky = 1:sampy*Ly
    ω  = 1:T

    lsgrid = GLMakie.labelslidergrid!(
        fig,
        ["kx", "ky", "ω"],
        [1:sampx*Lx, 1:sampy*Ly, 1:T]
    )
    fig[2, 1] = lsgrid.layout
    volslices = GLMakie.volumeslices!(ax, kx, ky, ω, sfactor)

    # See: http://makie.juliaplots.org/stable/plotting_functions/volumeslices.html
    sl_yz, sl_xz, sl_xy = lsgrid.sliders
    GLMakie.on(sl_yz.value) do v; volslices[:update_yz][](v) end
    GLMakie.on(sl_xz.value) do v; volslices[:update_xz][](v) end
    GLMakie.on(sl_xy.value) do v; volslices[:update_xy][](v) end

    GLMakie.set_close_to!(sl_yz, .5Lx)
    GLMakie.set_close_to!(sl_xz, .5Ly)
    GLMakie.set_close_to!(sl_xy, .5T)

    fig
end


text:
using Test, Lattices

@testset "Boundaries" begin
    @test Lattices.check_boundaries(Periodic(), Open()) == (Periodic(), Open())
    @test Lattices.check_boundaries(Helical(), Helical()) == (Helical(), Helical())

    @test Lattices.check_boundaries((Periodic(), Open())) == (Periodic(), Open())
    @test Lattices.check_boundaries((Helical(), Helical())) == (Helical(), Helical())

    @test_throws ArgumentError Lattices.check_boundaries(Periodic(), Helical())
end

@testset "Site IDs" begin
    # test 1D lattice
    for len in [1, 5, 10]
        l = Chain(len)
        for i in 1:len
            @test to_coordinate(l, to_site_id(l, Coordinate(i))) == Coordinate(i)
            @test to_site_id(l, to_coordinate(l, i)) == i
            @test to_coordinate(l, i) == Coordinate(i)
        end
    end

    # test a simple high-dimensional lattice
    l = HyperCubic((5,4,2,7))
    @test to_coordinate(l, to_site_id(l, Coordinate(1,1,1,1))) === Coordinate(1,1,1,1)  # origin
    @test to_coordinate(l, to_site_id(l, Coordinate(2,3,1,6))) === Coordinate(2,3,1,6)  # random bulk site
    @test to_coordinate(l, to_site_id(l, Coordinate(5,4,2,7))) === Coordinate(5,4,2,7)  # "last" site
    @test to_site_id(l, to_coordinate(l, 1)) == 1
    @test to_site_id(l, to_coordinate(l, nsites(l)÷10)) == nsites(l)÷10
    @test to_site_id(l, to_coordinate(l, nsites(l))) == nsites(l)
    @test to_coordinate(l, 1) == Coordinate(1,1,1,1)
    @test to_coordinate(l, nsites(l)) == Coordinate(5,4,2,7)

    @test_throws BoundsError to_coordinate(l, 0)
    @test_throws BoundsError to_coordinate(l, nsites(l)+1)
    @test_throws BoundsError to_site_id(l, Coordinate(0,0,0,0))
    @test_throws BoundsError to_site_id(l, Coordinate(2,3,4,6))
    @test_throws BoundsError to_site_id(l, Coordinate(5,4,2,8))
    @test_throws ArgumentError to_site_id(l, Coordinate(2,3,1,4,5))
    @test_throws ArgumentError to_site_id(l, Coordinate(5,4,2))

    # test case of non-trivial unit-cell
    l = Kagome((5,7))
    @test to_coordinate(l, to_site_id(l, Coordinate(1,1,1))) === Coordinate(1,1,1)  # origin
    @test to_coordinate(l, to_site_id(l, Coordinate(2,3,2))) === Coordinate(2,3,2)  # random bulk site
    @test to_coordinate(l, to_site_id(l, Coordinate(5,7,3))) === Coordinate(5,7,3)  # "last" site
    @test to_site_id(l, to_coordinate(l, 1)) == 1
    @test to_site_id(l, to_coordinate(l, nsites(l)÷6)) == nsites(l)÷6
    @test to_site_id(l, to_coordinate(l, nsites(l))) == nsites(l)
    @test to_coordinate(l, 1) == Coordinate(1,1,1)
    @test to_coordinate(l, nsites(l)) == Coordinate(5,7,3)

    @test_throws BoundsError to_coordinate(l, 0)
    @test_throws BoundsError to_coordinate(l, nsites(l)+1)
    @test_throws BoundsError to_site_id(l, Coordinate(0,0,0))
    @test_throws BoundsError to_site_id(l, Coordinate(2,3,4))
    @test_throws BoundsError to_site_id(l, Coordinate(5,7,4))
    @test_throws ArgumentError to_site_id(l, Coordinate(2,4,2,4))
    @test_throws ArgumentError to_site_id(l, Coordinate(2,4))
end

include("hypercubic.jl")
include("triangular.jl")
include("honeycomb.jl")
include("kagome.jl")


text:
# https://discourse.julialang.org/t/generate-random-value-from-a-given-function-out-of-box/5793/3
using Waterfall

# include("")
include(joinpath(WATERFALL_DIR,"src","includes.jl"))
include(joinpath(WATERFALL_DIR,"bin","io.jl"))

function define_permute(df, nperm::T; kwargs...) where T<:Integer
    N = size(df,1)-2
    lst = [[collect(1:N)]; random_permutation(1:N, min(factorial(N),nperm))]
    return define_permute(df, lst; kwargs...)
end

function define_permute(df::DataFrames.DataFrame, lst::T; fun::Function=identity, kwargs...) where T<:AbstractArray
    cascades = Dict()
    heights = Dict()

    for k in lst
        cascade = define_from(Cascade{Data}, df;
            permutation=k,
            nsample=100,
            correlate=true,
            permute=false,
            ncor=true,
            kwargs...,
        )
        push!(cascades, k => cascade)

        v1 = cumulative_v(cascade; shift=-1.0, kwargs...)
        v2 = cumulative_v(cascade; shift= 0.0, kwargs...)
        push!(heights, k => v2.-v1)
    end

    # Make dictionary of heights into one cascade.
    cascade = copy(cascades[first(lst)])
    data = collect_data(cascade)

    median = cat([collect_value(calculate(cascades[k], Statistics.quantile, 0.5)) for k in lst]...; dims=2)

    set_value!.(data, vectorize(median))
    return cascade
end


text:
module GPR

using StaticArrays
using LinearAlgebra
using Statistics
using Plots
using Optim
using ConstrainedDynamics
using ConstrainedDynamicsVis
using Rotations

export resetMechanism!
export updateMechanism!
export projectv!

export CState
export toState
export toStates
export setstates!

export MDCache
export MeanDynamics
export getμ

include("utils.jl")
include(joinpath("projections", "implicitProjection.jl"))
include("CState.jl")
include("mDynamics.jl")

end


text:
module StateVectorC

export _applyPauliX!, _applyHadamard!, _applyRX!, _applyRY!, _applyCX!
export applyPauliX, applyHadamard, applyRX, applyRY, applyCX

Base.@ccallable function apply1(gateLabel::String, sv::Ptr{ComplexF64}, num_elements::Int, wire::Int, inverse::Bool)::Cvoid
    sv_local = unsafe_wrap(Array, sv, num_elements; own=false)
    f = getfield(current_module(), gateLabel)
    f(sv_local, wire, inverse)
    ;
end
Base.@ccallable function apply1P(gateLabel::String, sv::Ptr{ComplexF64}, num_elements::Int, wire::Int, inverse::Bool, param::Float64)::Cvoid
    sv_local = unsafe_wrap(Array, sv, num_elements; own=false)
    f = getfield(current_module(), gateLabel)
    f(sv_local, wire, inverse, param)
    ;
end
Base.@ccallable function apply2(gateLabel::String, sv::Ptr{ComplexF64}, num_elements::Int, ctrl::Int, tgt::Int, inverse::Bool)::Cvoid
    sv_local = unsafe_wrap(Array, sv, num_elements; own=false)
    f = getfield(current_module(), gateLabel)
    f(sv_local, ctrl, tgt, inverse)
    ;
end
Base.@ccallable function apply2P(gateLabel::String, sv::Ptr{ComplexF64}, num_elements::Int, ctrl::Int, tgt::Int, inverse::Bool, param::Float64)::Cvoid
    sv_local = unsafe_wrap(Array, sv, num_elements; own=false)
    f = getfield(current_module(), gateLabel)
    f(sv_local, ctrl, tgt, inverse, param)
    ;
end

struct GateIndices
    internal::Vector{Int64}
    external::Vector{Int64}
end

mutable struct StateVector
    num_qubits::Int64
    data::Vector{ComplexF64}
end

function GateIndices(wires::Vector{Int64}, num_qubits::Int64)
    in = generateBitPatterns(wires, num_qubits)
    ex = generateBitPatterns(getIndicesAfterExclusion(wires, num_qubits), num_qubits)
    return GateIndices(in, ex)
end

function getIndicesAfterExclusion(indicesToExclude::Vector{Int64}, num_qubits::Int64)::Vector{Int64}
    indices = Vector{Int64}(1:1:num_qubits)
    filter!(x->x ∉ indicesToExclude, indices)
    return indices
end

function generateBitPatterns(qubitIndices::Vector{Int64}, num_qubits::Int64)::Vector{Int64}
    indices = zeros(1)
    for index_it = Iterators.reverse(qubitIndices)
        value = exp2(num_qubits - index_it)
        current_size = length(indices)
        for i = 1:1:current_size
            push!(indices, indices[i] + value)
        end
    end
    return indices
end

function _applyPauliX!(sv::StateVector,
                    indices::Vector{Int64},
                    externalIndices::Vector{Int64},
                    inverse::Bool = false)::Cvoid
    for ext_idx = externalIndices
        tmp1 = sv[1 + ext_idx + indices[1]]
        sv[1 + ext_idx + indices[1]] = sv[1 + ext_idx + indices[2]]
        sv[1 + ext_idx + indices[2]] = tmp
    end
    ;
end

function _applyHadamard!(sv::StateVector,
                        indices::Vector{Int64},
                        externalIndices::Vector{Int64},
                        inverse::Bool = false)::Cvoid
    for ext_idx = externalIndices
        v1 = sv[1 + ext_idx + indices[1]]
        v2 = sv[1 + ext_idx + indices[2]]

        sv[1 + ext_idx + indices[1]] = (1/sqrt(2))*(v1 + v2)
        sv[1 + ext_idx + indices[2]] = (1/sqrt(2))*(v1 - v2)
    end
    ;
end

function _applyRX!(sv::StateVector,
                    indices::Vector{Int64},
                    externalIndices::Vector{Int64},
                    param::Float64, inverse::Bool = false)::Cvoid
    c = cos(param / 2);
    js = (inverse == true) ?  -1im*sin(-param / 2) : 1im*sin(-param / 2);

    for ext_idx = externalIndices
        v1 = sv[1 + ext_idx + indices[1]]
        v2 = sv[1 + ext_idx + indices[2]]

        sv[1 + ext_idx + indices[1]] = c * v1 + js * v2
        sv[1 + ext_idx + indices[2]] = js * v1 + c * v2
    end
    ;
end

function _applyRY!( sv::StateVector,
                    indices::Vector{Int64},
                    externalIndices::Vector{Int64},
                    param::Float64, inverse::Bool = false)::Cvoid
    c = cos(param / 2);
    s = (inverse == true) ? -sin(param / 2) : sin(param / 2);

    for ext_idx = externalIndices
        v1 = sv[1 + ext_idx + indices[1]]
        v2 = sv[1 + ext_idx + indices[2]]

        sv[1 + ext_idx + indices[1]] = c * v1 - s * v2
        sv[1 + ext_idx + indices[2]] = s * v1 + c * v2
    end
    ;
end

function _applyCX!(sv::StateVector,
    indices::Vector{Int64},
    externalIndices::Vector{Int64},
    inverse::Bool = false)::Cvoid
    for ext_idx = externalIndices
        tmp1 = sv[1 + ext_idx + indices[3]]
        sv[1 + ext_idx + indices[3]] = sv[1 + ext_idx + indices[4]]
        sv[1 + ext_idx + indices[4]] = tmp
    end
    ;
end

Base.@ccallable function createSV(sv::StateVector, wire::Int64, inverse::Bool)::Cvoid
    gi = GateIndices([wire+1], length(sv));
    _applyPauliX!(sv, gi.internal, gi.external, inverse);
end

Base.@ccallable function applyPauliX(sv::StateVector, wire::Int64, inverse::Bool)::Cvoid
    gi = GateIndices([wire+1], length(sv));
    _applyPauliX!(sv, gi.internal, gi.external, inverse);
end

Base.@ccallable function applyHadamard(sv::StateVector, wire::Int64, inverse::Bool)::Cvoid
    gi = GateIndices([wire+1], length(sv));
    _applyHadamard!(sv, gi.internal, gi.external, inverse);
end

Base.@ccallable function applyRX(sv::StateVector, wire::Int64, inverse::Bool, param::Float64)::Cvoid
    gi = GateIndices([wire+1], length(sv));
    _applyRX!(sv, gi.internal, gi.external, param, inverse);
end

Base.@ccallable function applyRY(sv::StateVector, wire::Int64, inverse::Bool, param::Float64)::Cvoid
    gi = GateIndices([wire+1], length(sv));
    _applyRY!(sv, gi.internal, gi.external, param, inverse);
end

Base.@ccallable function applyCX(sv::StateVector, ctrl_wire::Int64, tgt_wire::Int64, inverse::Bool)::Cvoid
    gi = GateIndices([ctrl_wire+1, tgt_wire+1], length(sv));
    _applyCX!(sv, gi.internal, gi.external, inverse);
end

end

max_stars_repo_name: mlxd/PLext

max_issues_repo_name: mlxd/PLext

max_forks_repo_name: mlxd/PLext


text:
# Exact models

#src We generate a documentation page from this source file.
#src Let's tell the user:
#md # This file is automatically generated from source code.
#md # For usage examples refer to [Summary & Quick Examples](@ref exact_summary).

# ## Introduction
# The `ExactModel` is used to evaluate the objective exactly, without surrogate modelling
# (except for internal variable scaling).
# The derivatives are either user provided callbacks or can be deterimned using `ForwardDiff`
# or `FiniteDiff` automatically.

# ## Surrogate Model Interface Implementations

"""
    ExactModel( tfn, objf, diff_fn )

Exact Model type for evaluating the objective function `objf` directly.
Is instantiated by the corresponding `init_model` and `update_model` functions.
"""
struct ExactModel{
        M <: TransformerFn,
        O <: AbstractObjective,
        D <: DiffFn
    } <: SurrogateModel

    ## reference to a `TransformerFn` to have unscaling availabe:
    tfn :: M

    ## reference to objective(s) to evaluate
    objf :: O

    ## a `DiffFn` providing derivative information
    diff_fn :: Union{D,Nothing}
end

# The can determine the behavior of an `ExactModel` using `ExactConfig`:

"""
    ExactConfig(; gradients, jacobian = nothing, max_evals = typemax(Int64))

Configuration for an `ExactModel`.
`gradients` should be a vector of callbacks for the objective gradients **or**
a `Symbol`, either `:autodiff` or `fdm`, to define the differentiation method
to use on the objective.
Alternatively, a `jacobian` handle can be provided.
"""
@with_kw struct ExactConfig{
        G <: Union{Symbol, Nothing, AbstractVector{<:Function} },
        J <: Union{Nothing, Function}
    } <: SurrogateConfig

    gradients :: G = :autodiff
    ## alternative keyword, usage discouraged...
    jacobian :: J = nothing

    max_evals :: Int64 = typemax(Int64)

    @assert !(isnothing(gradients) && isnothing(jacobian)) "Provide either `gradients` or `jacobian`."
    @assert !(gradients isa Symbol) || gradients in [:autodiff, :fdm ] "`gradients` must be `:autodiff` or `:fdm`"
end

# There is no need for custom meta information:
struct ExactMeta <: SurrogateMeta end   # no construction meta data needed

# The remaining implementations are straightforward:
max_evals( emc :: ExactConfig ) = emc.max_evals
# We always deem the models fully linear:
fully_linear( em :: ExactModel ) = true
# They are not combinable to have individiual gradients availabe:
combinable( :: ExactConfig ) = false

# ## Construction

# When `cfg.gradients` is a `Symbol` we make use of the `DiffWrapper`s defined
# in `src/diff_wrappers.jl`:
function get_DiffFn( cfg :: ExactConfig{G,J}, objf :: AbstractObjective, tfn ) where{G<:Symbol,J}
    if cfg.gradients == :autodiff
        return AutoDiffWrapper( objf, tfn, nothing )
    elseif cfg.gradients == :fdm
        return FiniteDiffWrapper( objf, tfn, nothing );
    end
end

# Else we use a `GradWrapper`:
function get_DiffFn( cfg :: ExactConfig{G,J}, objf :: AbstractObjective, tfn) where{G,J}
    @assert length(cfg.gradients) == num_outputs(objf) "Provide as many gradient functions as the objective has outputs."
    return GradWrapper( tfn, cfg.gradients, cfg.jacobian )
end

# All "construction" work is done in the `_init_model` function:
function prepare_init_model(cfg ::ExactConfig, objf :: AbstractObjective,
    mop :: AbstractMOP, ::AbstractIterData, ::AbstractDB, :: AbstractConfig; kwargs...)
    return ExactMeta()
end

@doc "Return an ExactModel build from a VectorObjectiveFunction `objf`.
Model is the same inside and outside of criticality round."
function _init_model(cfg ::ExactConfig, objf :: AbstractObjective,
    mop :: AbstractMOP, ::AbstractIterData, ::AbstractDB, :: AbstractConfig, emeta :: ExactMeta; kwargs...)
    tfn = TransformerFn(mop)
    diff_fn = get_DiffFn( cfg, objf, tfn )
    em = ExactModel(tfn, objf, diff_fn )
    return em, emeta
end

# All the other functions simply return the input:
function update_model( em :: ExactModel, :: AbstractObjective, meta ::ExactMeta,
    ::AbstractMOP, :: AbstractIterData, ::AbstractDB, :: AbstractConfig;
    ensure_fully_linear :: Bool = false, kwargs... )
    return em, meta
end

function improve_model( em :: ExactModel, :: AbstractObjective, meta ::ExactMeta,
    ::AbstractMOP, :: AbstractIterData, ::AbstractDB, :: AbstractConfig;
    ensure_fully_linear :: Bool = false, kwargs ... )
    return em, meta
end

# ## Evaluation
@doc "Evaluate the ExactModel `em` at scaled site `x̂`."
function eval_models( em :: ExactModel, x̂ :: Vec )
    return eval_objf( em.objf, em.tfn, x̂ )
    ## using `eval_objf` will increase the evaluation count of `em.objf`
    ## That is why this count might be very high when using backtracking.
    ## eval_handle( em.objf )(x̂) would not increase the count.
end

@doc "Evaluate output `ℓ` of the ExactModel `em` at scaled site `x̂`."
function eval_models( em :: ExactModel, x̂ :: Vec, ℓ :: Int64)
    return eval_models(em,x̂)[ℓ]
end

@doc "Gradient vector of output `ℓ` of `em` at scaled site `x̂`."
function get_gradient( em :: ExactModel, x̂ :: Vec, ℓ :: Int64)
    return get_gradient( em.diff_fn, x̂, ℓ )
end

@doc "Jacobian Matrix of ExactModel `em` at scaled site `x̂`."
function get_jacobian( em :: ExactModel, x̂ :: Vec )
    return get_jacobian( em.diff_fn, x̂ );
end

# ## [Summary & Quick Examples](@id exact_summary)

# ```julia
# using Morbit
# Morbit.print_all_logs()
# mop = MixedMOP(3)
#
# f1 = x -> sum( ( x .- 1 ).^2
# f2 = x -> sum( ( x .+ 1 ).^2
# g1 = x -> 2 .* ( x .- 1 )
# g2 = x -> 2 .* ( x .+ 1 )
#
# add_objective!( mop, f1, ExactConfig(; gradients = [g1,]) )
# add_objective!( mop, f2, ExactConfig(; gradients = [g2,]) )
#
# x_fin, f_fin, _ = optimize( mop, [-π, ℯ, 0])
# ```


text:
using Base.depwarn

typealias EmptyStreamSource EmptyStream

Base.@deprecate EmptyStreamSource EmptyStream
Base.@deprecate seekforward(stream::BufferedInputStream, n::Integer) skip(stream, n)


text:
"""
$(TYPEDEF)

Hypograph of perspective function of sum-log cone of dimension `dim`.

    $(FUNCTIONNAME){T}(dim::Int, use_dual::Bool = false)
"""
mutable struct HypoPerLog{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    ϕ::T
    ζ::T
    tempw::Vector{T}

    function HypoPerLog{T}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real}
        @assert dim >= 3
        cone = new{T}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        return cone
    end
end

function setup_extra_data!(cone::HypoPerLog{T}) where {T <: Real}
    d = cone.dim - 2
    cone.tempw = zeros(T, d)
    return cone
end

get_nu(cone::HypoPerLog) = cone.dim

function set_initial_point!(arr::AbstractVector, cone::HypoPerLog)
    (arr[1], arr[2], w) = get_central_ray_hypoperlog(cone.dim - 2)
    @views arr[3:end] .= w
    return arr
end

function update_feas(cone::HypoPerLog{T}) where {T <: Real}
    @assert !cone.feas_updated
    v = cone.point[2]
    @views w = cone.point[3:end]

    if v > eps(T) && all(>(eps(T)), w)
        u = cone.point[1]
        cone.ϕ = sum(log(wi / v) for wi in w)
        cone.ζ = v * cone.ϕ - u
        cone.is_feas = (cone.ζ > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::HypoPerLog{T}) where {T <: Real}
    u = cone.dual_point[1]
    @views w = cone.dual_point[3:end]

    if all(>(eps(T)), w) && u < -eps(T)
        v = cone.dual_point[2]
        sumlog = sum(log(w_i / -u) for w_i in w)
        return (v - u * (sumlog + length(w)) > eps(T))
    end

    return false
end

function update_grad(cone::HypoPerLog)
    @assert cone.is_feas
    u = cone.point[1]
    v = cone.point[2]
    @views w = cone.point[3:end]
    d = length(w)
    ζ = cone.ζ
    g = cone.grad

    g[1] = inv(cone.ζ)
    g[2] = -(cone.ϕ - d) / ζ - inv(v)
    vζi1 = -1 - v / ζ
    @. g[3:end] = vζi1 / w

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::HypoPerLog)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    v = cone.point[2]
    @views w = cone.point[3:end]
    H = cone.hess.data
    ζ = cone.ζ
    d = length(w)
    σζi = (cone.ϕ - d) / ζ
    vζi = v / ζ

    # u, v
    H[1, 1] = ζ^-2
    H[1, 2] = -σζi / ζ
    H[2, 2] = v^-2 + abs2(σζi) + d / (ζ * v)

    # u, v, w
    vζi2 = -vζi / ζ
    c1 = ((cone.ϕ - d) * vζi - 1) / ζ
    @. H[1, 3:end] = vζi2 / w
    @. H[2, 3:end] = c1 / w

    # w, w
    @inbounds for j in eachindex(w)
        j2 = 2 + j
        wj = w[j]
        c2 = vζi / wj
        for i in 1:(j - 1)
            H[2 + i, j2] = vζi / w[i] * c2
        end
        H[j2, j2] = (c2 + c2 * vζi + inv(wj)) / wj
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::HypoPerLog,
    )
    v = cone.point[2]
    @views w = cone.point[3:end]
    d = length(w)
    ζ = cone.ζ
    σ = cone.ϕ - d
    vζi1 = (v + ζ) / ζ
    rwi = cone.tempw

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @. @views rwi = arr[3:end, j] / w

        qζi = q / ζ
        c1 = sum(rwi) / ζ
        # ∇ϕ[r] = v * c1
        c2 = (v * c1 - p / ζ + σ * qζi) / ζ
        c3 = c2 * v - qζi

        prod[1, j] = -c2
        prod[2, j] = c2 * σ - c1 + (qζi * d + q / v) / v
        @. prod[3:end, j] = (c3 + vζi1 * rwi) / w
    end

    return prod
end

function update_inv_hess(cone::HypoPerLog)
    @assert cone.grad_updated
    @assert !cone.inv_hess_updated
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    v = cone.point[2]
    @views w = cone.point[3:end]
    d = length(w)
    Hi = cone.inv_hess.data
    ζ = cone.ζ
    ϕ = cone.ϕ
    ζv = ζ + v
    ζζvi = ζ / ζv
    c1 = v / (ζv + d * v) * v
    c2 = c1 / ζv

    # u, v
    Hi12 = Hi[1, 2] = c1 * (ζv * ϕ - d * ζ)
    Hi[1, 1] = ζζvi * d * v^2 + ζ^2 + (ϕ - ζζvi * d) * Hi12
    Hi[2, 2] = c1 * ζv

    # u, v, w
    c3 = (v * ζ + Hi12) / ζv
    @. Hi[1, 3:end] = c3 * w
    @. Hi[2, 3:end] = c1 * w

    # w, w
    @inbounds for j in eachindex(w)
        j2 = 2 + j
        wj = w[j]
        c4 = c2 * wj
        for i in 1:j
            Hi[2 + i, j2] = c4 * w[i]
        end
        Hi[j2, j2] += wj * ζζvi * wj
    end

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::HypoPerLog,
    )
    @assert cone.grad_updated
    v = cone.point[2]
    @views w = cone.point[3:end]
    d = length(w)
    ζ = cone.ζ
    ϕ = cone.ϕ
    ζv = ζ + v
    ζζvi = ζ / ζv
    c1 = ζv / v + d
    rw = cone.tempw

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @. @views rw = arr[3:end, j] * w

        trrw = sum(rw)
        c2 = v * (ζv * (ϕ * p + q) - d * ζ * p + trrw) / c1
        c3 = (c2 + ζ * v * p) / ζv

        prod[1, j] = ζ * ((v * (d * p * v + trrw) - d * c2) / ζv + ζ * p) + ϕ * c2
        prod[2, j] = c2
        @. prod[3:end, j] = (c3 + ζζvi * rw) * w
    end

    return prod
end

function dder3(cone::HypoPerLog{T}, dir::AbstractVector{T}) where {T <: Real}
    @assert cone.grad_updated
    v = cone.point[2]
    @views w = cone.point[3:end]
    dder3 = cone.dder3
    p = dir[1]
    q = dir[2]
    ζ = cone.ζ
    d = length(w)
    σ = cone.ϕ - d
    viq = q / v
    vζi1 = (v + ζ) / ζ
    rwi = cone.tempw

    @. @views rwi = dir[3:end] / w
    tr1 = sum(rwi)
    tr2 = sum(abs2, rwi)

    χ = (-p + σ * q + tr1 * v) / ζ
    c1 = (viq * (2 * tr1 - viq * d) - tr2) / (2 * ζ)
    c2 = (abs2(χ) - v * c1) / ζ
    c3 = -q * χ / ζ
    c4 = (χ * v - q) / ζ
    c5 = c3 + c2 * v

    dder3[1] = -c2
    dder3[2] = c2 * σ + (abs2(viq) - d * c3 - c4 * tr1) / v - tr2 / ζ - c1
    @. dder3[3:end] = (c5 + rwi * (c4 + vζi1 * rwi)) / w

    return dder3
end

# see analysis in
# https://github.com/lkapelevich/HypatiaSupplements.jl/tree/master/centralpoints
function get_central_ray_hypoperlog(d::Int)
    if d <= 10
        # lookup points where x = f'(x)
        return central_rays_hypoperlog[d, :]
    end
    # use nonlinear fit for higher dimensions
    x = inv(d)
    if d <= 70
        u = 4.657876 * x ^ 2 - 3.116192 * x + 0.000647
        v = 0.424682 * x + 0.553392
        w = 0.760412 * x + 1.001795
    else
        u = -3.011166 * x - 0.000122
        v = 0.395308 * x + 0.553955
        w = 0.837545 * x + 1.000024
    end
    return [u, v, w]
end

const central_rays_hypoperlog = [
    -0.827838387  0.805102007  1.290927686
    -0.689607388  0.724605082  1.224617936
    -0.584372665  0.68128058  1.182421942
    -0.503499342  0.65448622  1.153053152
    -0.440285893  0.636444224  1.131466926
    -0.389979809  0.623569352  1.114979519
    -0.349255921  0.613978276  1.102013921
    -0.315769104  0.606589839  1.091577908
    -0.287837744  0.600745284  1.083013
    -0.264242734  0.596019009  1.075868782
    ]


text:
using Test
using DataFrames, TableOperations, Tables, Random
using Recommenders:
    Movielens100k,
    load_dataset,
    ratio_split,
    ItemkNN,
    evaluate_u2i,
    MeanPrecision,
    MeanRecall,
    MeanNDCG

ml100k = Movielens100k()
download(ml100k)
rating, _, _ = load_dataset(ml100k)

Random.seed!(1234);
train_table, test_table = ratio_split(rating, 0.8)

prec10 = MeanPrecision(10)
recall10 = MeanRecall(10)
ndcg10 = MeanNDCG(10)
metrics = [prec10, recall10, ndcg10]

model = ItemkNN(10, 0.001, :bm25, true, true, true)
result = evaluate_u2i(
    model,
    train_table,
    test_table,
    metrics,
    10,
    col_user = :userid,
    col_item = :movieid,
    col_rating = :rating,
    drop_history = true,
)

@test !(model.similarity === nothing)

@test result[:ndcg10] >= 0
@test result[:precision10] >= 0
@test result[:recall10] >= 0

pred_u2i = predict_u2i(model, 1, 10)
@test length(pred_u2i) == 10

pred_u2i = predict_u2i(model, [1, 10], 10)
for p in pred_u2i
    @test length(p) == 10
end

pred_i2i = predict_i2i(model, 1, 10)
@test length(pred_i2i) == 10

pred_i2i = predict_i2i(model, [1, 10], 10)
for p in pred_i2i
    @test length(p) == 10
end

tempfile = tempname()
save_model(model, tempfile)
@test isfile(tempfile)
model = load_model(tempfile)
@test model.similarity !== nothing


text:
function test_projections()
    Dists = [Uniform(0,1), Normal(0,1)]
    C1s   = [0.001, 0.1, 1, 10, 100]
    C2s   = [0.001, 0.1, 1, 10, 100]
    Ks    = 1:10
    δ0s   = vcat(0:0.1:1, 10, 100)

    n = rand(100:200)
    m = rand(50:150)

    @testset "projections" begin
        @testset "α0 ∼ $(nameof(typeof(distα))), β0 ∼ $(nameof(typeof(distβ)))" for (distα, distβ) in Iterators.product(Dists, Dists)
            α0 = rand(distα, n)
            β0 = rand(distα, m)

            @testset "δ0 = $δ0, C1 = $C1, C2 = $C2" for (δ0, C1, C2) in Iterators.product(δ0s, C1s, C2s)
                test_projection(α0, β0, δ0, C1, C2)
            end
            @testset "δ0 = $δ0" for δ0 in δ0s
                test_projection(α0, β0, δ0)
            end
            @testset "C = $C, K = $K" for (C, K) in Iterators.product(C1s, Ks)
                test_projection(α0, β0, C, K)
            end
            @testset "K = $K" for K in Ks
                test_projection(α0, β0, K)
            end
        end
    end
end

function objective(α::AbstractVector, β::AbstractVector, δ::Real, α0::AbstractVector, β0::AbstractVector, δ0::Real)
    sum(abs2, α - α0)/2 + sum(abs2, β - β0)/2 + sum(abs2, δ - δ0)/2
end

function objective(α::AbstractVector, β::AbstractVector, α0::AbstractVector, β0::AbstractVector)
    sum(abs2, α - α0)/2 + sum(abs2, β - β0)/2
end

function test_projection(α0::AbstractVector, β0::AbstractVector, δ0::Real, C1::Real, C2::Real; atol::Real = 1e-5)
    α1, β1, δ1 = AccuracyAtTopKernels.projection(α0, β0, δ0, C1, C2)
    α2, β2, δ2 = AccuracyAtTopKernels.projection_exact(α0, β0, δ0, C1, C2)

    @testset "patmat hinge" begin
        @testset "feasibility" begin
            @test abs(sum(α1) - sum(β1)) <= atol
            @test minimum(α1) >= - atol
            @test maximum(α1) <= C1 + atol
            @test minimum(β1) >= - atol
            @test maximum(β1) <= C2*δ1 + atol
        end
        @testset "optimality" begin
            @test objective(α1, β1, δ1, α0, β0, δ0) <= objective(α2, β2, δ2, α0, β0, δ0) + atol
        end
    end
end

function test_projection(α0::AbstractVector, β0::AbstractVector, δ0::Real; atol::Real = 1e-5)
    α1, β1, δ1 = AccuracyAtTopKernels.projection(α0, β0, δ0)
    α2, β2, δ2 = AccuracyAtTopKernels.projection_exact(α0, β0, δ0)

    @testset "patmat quadratic" begin
        @testset "feasibility" begin
            @test abs(sum(α1) - sum(β1)) <= atol
            @test minimum(α1) >= - atol
            @test minimum(β1) >= - atol
            @test minimum(δ1) >= - atol
        end
        @testset "optimality" begin
            @test objective(α1, β1, δ1, α0, β0, δ0) <= objective(α2, β2, δ2, α0, β0, δ0) + atol
        end
    end
end

function test_projection(α0::AbstractVector, β0::AbstractVector, C::Real, K::Integer; atol::Real = 1e-5)
    α1, β1 = AccuracyAtTopKernels.projection(α0, β0, C, K)
    α2, β2 = AccuracyAtTopKernels.projection_exact(α0, β0, C, K)

    @testset "toppushk hinge" begin
        @testset "feasibility" begin
            @test abs(sum(α1) - sum(β1)) <= atol
            @test minimum(α1) >= - atol
            @test maximum(α1) <= C + atol
            @test minimum(β1) >= - atol
            @test maximum(β1) <= sum(α1)/K + atol
        end
        @testset "optimality" begin
            @test objective(α1, β1, α0, β0) <= objective(α2, β2, α0, β0) + atol
        end
    end
end

function test_projection(α0::AbstractVector, β0::AbstractVector, K::Integer; atol::Real = 1e-5)
    α1, β1 = AccuracyAtTopKernels.projection(α0, β0, K)
    α2, β2 = AccuracyAtTopKernels.projection_exact(α0, β0, K)

    @testset "toppushk quadratic" begin
        @testset "feasibility" begin
            @test abs(sum(α1) - sum(β1)) <= atol
            @test minimum(α1) >= - atol
            @test minimum(β1) >= - atol
            @test maximum(β1) <= sum(α1)/K + atol
        end
        @testset "optimality" begin
            @test objective(α1, β1, α0, β0) <= objective(α2, β2, α0, β0) + atol
        end
    end
end


text:
@testset "a† and a" begin
    dim = DIM

    @test create!(VacuumState(dim=dim)) ≈ SinglePhotonState(dim=dim)
    @test annihilate!(SinglePhotonState(dim=dim)) ≈ VacuumState(dim=dim)
    @test create!(VacuumState(dim=dim, rep=StateMatrix)) ≈ SinglePhotonState(dim=dim, rep=StateMatrix)
    @test annihilate!(SinglePhotonState(dim=dim, rep=StateMatrix)) ≈ VacuumState(dim=dim, rep=StateMatrix)

    @test create(VacuumState(dim=dim)) ≈ SinglePhotonState(dim=dim)
    @test annihilate(SinglePhotonState(dim=dim)) ≈ VacuumState(dim=dim)
    @test create(VacuumState(dim=dim, rep=StateMatrix)) ≈ SinglePhotonState(dim=dim, rep=StateMatrix)
    @test annihilate(SinglePhotonState(dim=dim, rep=StateMatrix)) ≈ VacuumState(dim=dim, rep=StateMatrix)
end

@testset "α and ξ" begin
    @test repr(ComplexVec(2., π/4)) == "ComplexVec{Float64}(2.0exp(-$(π/4)im))"
    @test QSB.z(α(2, π/4)) ≈ 2 * exp(-im * π/4)
    @test QSB.z(ξ(2, π/4)) ≈ 2 * exp(-im * π/4)
end

@testset "Displacement" begin
    dim = DIM
    r = 2.
    θ = π/4

    @test displace!(VacuumState(dim=dim), α(r, θ)).v ≈ exp(
        QSB.z(α(r, θ)) * Creation(dim=dim) -
        QSB.z(α(r, θ))' * Annihilation(dim=dim)
    ) * VacuumState().v
    @test displace!(VacuumState(dim=dim, rep=StateMatrix), α(r, θ)).𝛒 ≈ exp(
        QSB.z(α(r, θ)) * Creation(dim=dim) -
        QSB.z(α(r, θ))' * Annihilation(dim=dim)
    ) * VacuumState(rep=StateMatrix).𝛒 * exp(
        QSB.z(α(r, θ)) * Creation(dim=dim) -
        QSB.z(α(r, θ))' * Annihilation(dim=dim)
    )'
end

@testset "squeezing" begin
    dim = DIM
    r = 2.
    θ = π/4

    @test squeeze!(VacuumState(dim=dim), α(r, θ)).v ≈ exp(
        0.5 * QSB.z(ξ(r, θ))' * Annihilation(dim=dim)^2 -
        0.5 * QSB.z(ξ(r, θ)) * Creation(dim=dim)^2
    ) * VacuumState().v
    @test squeeze!(VacuumState(dim=dim, rep=StateMatrix), α(r, θ)).𝛒 ≈ exp(
        0.5 * QSB.z(ξ(r, θ))' * Annihilation(dim=dim)^2 -
        0.5 * QSB.z(ξ(r, θ)) * Creation(dim=dim)^2
    ) * VacuumState(rep=StateMatrix).𝛒 * exp(
        0.5 * QSB.z(ξ(r, θ))' * Annihilation(dim=dim)^2 -
        0.5 * QSB.z(ξ(r, θ)) * Creation(dim=dim)^2
    )'
end

@testset "measurement" begin
    ψₙs = QSB.ψₙ.(big(0):big(DIM-1), 2., 3.)
    @test QSB.𝛑̂(2, 3, dim=DIM) ≈ ψₙs * ψₙs'
end

@testset "Gaussian state" begin
    𝐚 = rand(10, 10)
    𝐛 = rand(10, 10)

    @test QSB.tr_mul(𝐚, 𝐛) ≈ tr(𝐚 * 𝐛)

    state = SqueezedThermalState(ξ(1., π/4), 0.5)

    @test QSB.create_μ(state) ≈ tr(Creation(dim=state.dim) * state.𝛒)
    @test QSB.create²_μ(state) ≈ tr(Creation(dim=state.dim)^2 * state.𝛒)
    @test QSB.annihilate_μ(state) ≈ tr(Annihilation(dim=state.dim) * state.𝛒)
    @test QSB.annihilate²_μ(state) ≈ tr(Annihilation(dim=state.dim)^2 * state.𝛒)
    @test QSB.create_annihilate_μ(state) ≈ tr(Creation(dim=state.dim) * Annihilation(dim=state.dim) * state.𝛒)
end


text:
using PyRhodium
using Roots
using Distributions
using DataFrames
using Random

function lake_problem(;pollution_limit = nothing,
         b = 0.42,       # decay rate for P in lake (0.42 = irreversible)
         q = 2.0,        # recycling exponent
         μ = 0.02,       # mean of natural inflows
         σ = 0.001,      # standard deviation of natural inflows
         α = 0.4,        # utility from pollution
         δ = 0.98,       # future utility discount rate
         nsamples = 100) # monte carlo sampling of natural inflows)

    Pcrit = fzero(x->x^q / (1 + x^q) - b * x, 0.01, 1.5)
    nvars = length(pollution_limit)
    X = zeros(nvars)
    average_daily_P = zeros(nvars)
    reliability = 0.0

    d = LogNormal(log(μ^2 / sqrt(σ^2 + μ^2)), sqrt(log(1.0 + σ^2 / μ^2)))

    natural_inflows = zeros(nvars)

    for i in 1:nsamples
        X[1] = 0.0

        rand!(d, natural_inflows)

        for t in 2:nvars
            X[t] = (1 - b) * X[t - 1] + X[t - 1]^q / (1 + X[t - 1]^q) + pollution_limit[t] + natural_inflows[t]
            average_daily_P[t] += X[t] / nsamples
        end

        reliability += sum(X .< Pcrit) / (nsamples * nvars)
    end

    max_P = maximum(average_daily_P)
    utility = sum(α .* pollution_limit .* δ.^collect(1:nvars))
    intertia = sum(diff(pollution_limit) .> -0.02) / (nvars - 1)

    return max_P, utility, intertia, reliability
end

m = Model(lake_problem)

set_parameters!(m, [Parameter("pollution_limit"),
                    Parameter("b"),
                    Parameter("q"),
                    Parameter("mean"),
                    Parameter("stdev"),
                    Parameter("delta")])

set_responses!(m, [Response("max_P", :MINIMIZE),
                   Response("utility", :MAXIMIZE),
                   Response("inertia", :MAXIMIZE),
                   Response("reliability", :MAXIMIZE)])

set_levers!(m, [RealLever("pollution_limit", 0.0, 0.1, length = 100)])

# output = optimize(m, "NSGAII", 1000)
output = optimize(m, "NSGAII", 100)

println("Found $(length(output)) optimal policies!")

# fig = scatter2d(output)

policy = output[5]
policy = findmax(output, :reliability)
policies = find(output, "utility > 0.5");

policy = named_tuple(findmax(output, :reliability))

println("Max Phosphorus in Lake: ", policy.max_P)
println("Utility:                ", policy.utility)
println("Inertia:                ", policy.inertia)
println("Reliability:            ", policy.reliability)

df = DataFrame(output)
arr = collect(output);


text:
#
# Returns the index in the natural_aminoacids list of the amino acid name given
#
function which_natural_aminoacid(atom::Atom)

    # If the residue name doesn't have at least three letters, this is not a protein atom
    if length(atom.resname) < 3
        return 0
    end

    # To take into account alternate conformations, such as "AGLY", "BGLY", etc.
    l = length(atom.resname)
    name = atom.resname[l-2:l]

    # Return the index of this amino acid in the amino acid list, or nothing
    i = findfirst(aa -> aa.three_letter_code == name, natural_aminoacids)
    if i == nothing
        return 0
    end

    return i
end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
module Datasets

using DrWatson
using UCI
using MLDatasets
using DataDeps
using StatsBase
using DelimitedFiles
using Random
using CSV
using DataFrames
using Flux # for one-hot encoding
using NPZ
using Images

export load_data

include("datadeps_init.jl")
include("tabular.jl")
include("img.jl")
include("basics.jl")

end


text:
"""
	aggregate2

Aggregate data to yearly, monthly, daily, hourly or ONE minute samples.
Simple aggregation = sum, mean, maximum or minimum can be applied during re-sampling.
Missing (NaN) data can be either kept or replaced during re-sampling.

**Input**
* data: DataFrame where at least one column contains DateTime
* timecol: column containing DateTime. Default column name = :datetime
* resol: output (date)time resolution, e.g. Dates.Hour(1) (default).
> Dates.Hour(2) is not allowed (only ONE hour|minute|...)

* fce: function to by applied = sum (default)| mean | maximum | minimum
> To remove NaNs use `x->sum(filter(!isnan,x))` function

**Output**
* resampled dataframe with all input columns
> **Warning** the output dataframe can be re-arranged!

**Example**
```
data = DataFrame(Temp=[10,11,12,14],
       datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
         DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
dataa = aggregate2(data,resol=Dates.Day(1),fce=x->minimum(collect(skipmissing(x))));
```
"""
function aggregate2(data::DataFrame;
					timecol=:datetime,resol=Dates.Hour(1),
					fce=sum,rename=false)
	dfc,datestringcol = addTimeString(data,resol,timecol);
	# do not use datetime for aggregation (not supported)
	useonly = allexcept(names(dfc),timecol);
	dfc = aggregate(dfc[useonly],:datestringcol,fce)
	# Add back back datetime and remove date-string used for aggregation
	dfc[timecol] = DateTime.(dfc[:datestringcol],datestringcol);
	# make suree that type is DateTime
	if eltype(dfc[timecol]) != DateTime
		dfc[timecol] = convert(Vector{DateTime},dfc[timecol]);
	end
	delete!(dfc,:datestringcol);
	if rename==true
		for i in names(dfc)
			temp = string(i) |> x-> findlast("_",x);
			if typeof(temp)!=Nothing
				rename!(dfc,i=>Symbol(string(i)[1:collect(temp)[1]-1]));
			end
		end
	end
	return dfc
end

"""
	time2pattern(resol)

Convert DateTime hour, minute,... to string pattern, e.g. "yyyymmhh" or "yyyymmddhhmm"

**Input**
* resol: output (date)time resolution, e.g. Dates.Hour(1). (Dates.Hour(2) is not allowed)

**Output**
* datetime string

**Example**
```
s = time2pattern(Dates.Hour(1));
```
"""
function time2pattern(resol)
	if resol==Dates.Second(1)
		return "yyyymmddHHMMSS"
	elseif resol==Dates.Minute(1)
		return "yyyymmddHHMM"
	elseif resol==Dates.Hour(1)
		return "yyyymmddHH"
	elseif resol==Dates.Day(1)
		return "yyyymmdd"
	elseif resol=="month";#Dates.Month(1)
		return "yyyymm"
	elseif resol=="year"#Dates.Year(1)
		return "yyyy"
	end
end

"""
	time2regular(data,timecol,resol)

Resample data to regular time sampling (missing filled with NA values)

**Input**
* data: DataFrame where at least one column contains DateTime
* timecol: column containing DateTime. Default column name = :datetime
* resol: output (date)time resolution, e.g. Dates.Hour(1)

**Output**
* resampled dataframe with all input columns

**Example**
```
data = DataFrame(Temp=[10,11,12,14],
			datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
			DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
reg_sample = time2regular(data,timecol=:datetime,resol=Dates.Hour(1))

```
"""
function time2regular(data::DataFrame;timecol=:datetime,resol=Dates.Hour(1))
	# create dataframe with regular sampling
	dfr = DataFrame(datetime=collect(data[timecol][1]:resol:data[timecol][end]));
	# joint original and regular sampled dataframe
	reg_sample = join(data,dfr,on=timecol,kind=:right);
	# Sort accoring to time
	return sort!(reg_sample, timecol);
end

"""
	isregular(timevec)

Check if the data/time vector is regularly sampled

**Input**
`timevec`: Vector with DateTime

**Output**
`true` for regularly sampled data

**Example**
```
timevec = [DateTime(2010,1,1,1,0,0),
				 DateTime(2010,1,1,2,0,0),
				 DateTime(2010,1,1,3,0,0),
				 DateTime(2010,1,1,5,0,0),# fourth hour is missing =>not regular
				 DateTime(2010,1,1,6,0,0)];
out = isregular(timevec); # will return false
```
"""
function isregular(timevec::Vector{DateTime})
	timediff = diff(Dates.value.(timevec));
	if all(timediff[1] .== timediff)
		return true;
	else
		return false;
	end
end

"""
	getresolution(timevec;timecol)
Get time resolution in `Dates` format
NOT tested yet.

**Input:**
`timevec`: vector with DateTime values

**Output:**
time resolution in `Dates` format. Will return one of following:
Dates.Millisecond(X) | Dates.Second(X) | Dates.Minute(X) | Dates.Hour(X) |
Dates.Day(X),  where X <0,10>

**Example:**
```
timevec = [DateTime(2010,1,1,1,0,0),
				 DateTime(2010,1,1,2,0,0),
				 DateTime(2010,1,1,3,0,0)];
resol = getresolution(timevec); # should return Dates.Hour(1);
```
"""
function getresolution(timevec::Vector{DateTime})
	dt = timevec[2]-timevec[1];
	if dt < Dates.Millisecond(1000) # millisecond
		return dt;
	elseif dt < Dates.Millisecond(60*1000) # seconds
		return Dates.Second(dt)
	elseif dt < Dates.Millisecond(60*60*1000) # minute
		return Dates.Minute(dt)
	elseif dt < Dates.Millisecond(24*60*60*1000) # hour
		return Dates.Hour(dt)
	elseif dt < Dates.Millisecond(24*60*60*60*1000) # day
		return Dates.Day(1)
	else
		return dt;
	end
end

"""
	allexept(name)

Auxiliary function to return all column numbers except for input name

**Input**
* head: names collection (=names(dataframe))
* name: name/key to be excluded

**Output**
* list of indices
"""
function allexcept(head,name)
	index = Vector{Int64}();
	for (i,val) in enumerate(head)
		if val != name
			push!(index,i)
		end
	end
	return index;
end

"""
Auxiliary function to prepare dataframe for aggregate2 and timejoin
"""
function addTimeString(data,resol,timecol)
	# Create a copy of the DateFrame for manipulation
	dfc = deepcopy(data);
	# Convert input resolution to string pattern + apply
	datestringcol = time2pattern(resol);
	# Append column with string date/time
	dfc[:datestringcol] = Dates.format.(data[timecol],datestringcol);
	return dfc,datestringcol
end

"""
	cut2interval!(datain,timestart,timeend;)
Cut input dataframe to set time interval

**Input**
* datain: input dataframe containing :datetime (DateTime)
* timestart: starting time (DateTime)
* timeend: end time (DateTime)
* keepedges: either keep (true) or remove (false) the time epoch at (timestart,timeend)

> nothing will be deleted if start and end point are outside of input time

**Example**
```
datatest = DataFrame(datetime=collect(DateTime(2000,1,1):Dates.Day(1):DateTime(2000,1,12)),
					 somevalues=collect(1:1:12))
timestart,timeend = DateTime(2000,1,2),DateTime(2000,1,11);
cut2interval!(datatest,starttime,endtime,keepedges=(false,true))
```
"""
function cut2interval!(datain::DataFrame,timestart::DateTime,timeend::DateTime;
						keepedges::Tuple{Bool,Bool}=(true,true))
	rf = ResampleAndFit.findinterval(datain[:datetime],timestart,timeend,
									includetime=true);
	if length(rf) > 0 # at least 1 point must be found to delete rows
		r1 = keepedges[1] ? collect(1:rf[1]-1) : collect(1:rf[1]);
		r2 = keepedges[2] ? collect(rf[end]+1:size(datain,1)) : collect(rf[end]:size(datain,1));
		if length(r1) > 0 || length(r2) > 0
			deleterows!(datain,vcat(r1,r2));
		end
	end
end


text:
#=
# 3D Poisson, Dirichlet bc
# Uses Dendro imported as a custom gen target
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("poissondendro");

useLog("poissondendrolog", level=3)

# default values (max_depth=6, wavelet_tol = 0.1, partition_tol = 0.3, solve_tol = 1e-6, max_iters = 100)
generateFor("target_dendro_cg.jl", params=(7, 0.01, 0.3, 0.000001, 1000))

domain(3)
functionSpace(order=4)

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, "0")

f = coefficient("f", "-30*pi*pi*sin(5*pi*x)*sin(2*pi*y)*sin(pi*z)")
weakForm(u, "-dot(grad(u),grad(v)) - f*v")

build_octree_with(f);
solve(u);

finalize_femshop()

max_stars_repo_name: paralab/femshop

max_issues_repo_name: paralab/femshop

max_forks_repo_name: paralab/femshop


text:
module LossFunctions
#-----------------------------------------------------------
#PURPOSE
#	Colin T. Bowers module for loss functions
#NOTES
#	To add a new loss function, provide a new type definition and any constructors under the TYPE DEFINTIIONS heading, then add an appropriate loss method for scalar inputs.
#	If you want to aggregate your loss series in some way, you should probably look at using the Distances package before thinking about this package
#LICENSE
#	MIT License (see github repository for more detail: https://github.com/colintbowers/LossFunctions.jl.git)
#-----------------------------------------------------------

#Load any entire modules that are needed (use import ModuleName1, ModuleName2, etc)

#Load any specific variables/functions that are needed (use import ModuleName1.FunctionName1, ModuleName2.FunctionName2, etc)
import 	Base.string,
		Base.show

#Specify the variables/functions to export (use export FunctionName1, FunctionName2, etc)
export 	LossFunction, #Abstract supertype
		SymmetricLossFunction, #Abstract subtypes
		AsymmetricLossFunction, #Abstract subtypes
		SquaredLoss,
		AbsoluteLoss,
		MinkowskiLoss,
		HuberLoss,
		SquaredLogLoss,
		AbsoluteLogLoss,
		SquaredPropLoss,
		AbsolutePropLoss,
		QLIKE,
		QLIKEReverse,
		loss, #Generic method for calculating loss, ie L(x, y), or L(e) for the case e = x - y
		lossdiff #Generic method for calculating a loss differential, ie L(x-BaseCase) - L(y-BaseCase)

#******************************************************************************

#----------------------------------------------------------
#SET CONSTANTS FOR MODULE
#----------------------------------------------------------
const singleInputLossMethods = Union{SquaredLoss, AbsoluteLoss, MinkowskiLoss, HuberLoss}::Union

#----------------------------------------------------------
#LOSS FUNCTION TYPES
#----------------------------------------------------------
#Abstract types
abstract LossFunction
abstract SymmetricLossFunction <: LossFunction
abstract AsymmetricLossFunction <: LossFunction
#Loss function types
type DummyLoss <: LossFunction ; end #dummy type that does nothing
type SquaredLoss <: SymmetricLossFunction ; end
type AbsoluteLoss <: SymmetricLossFunction ; end
type MinkowskiLoss <: SymmetricLossFunction
	p::Float64
	MinkowskiLoss(p::Float64) = (p <= 0) ? error("Minkowski parameter must be strictly positive") : new(p)
end
MinkowskiLoss() = MinkowskiLoss(2.0) #Default is SquaredLoss
MinkowskiLoss(p::Number) = MinkowskiLoss(convert(Float64, p))
type SquaredLogLoss <: SymmetricLossFunction ; end
type AbsoluteLogLoss <: SymmetricLossFunction ; end
type SquaredPropLoss <: AsymmetricLossFunction ; end
type AbsolutePropLoss <: AsymmetricLossFunction ; end
type QLIKE <: AsymmetricLossFunction ; end
type QLIKEReverse <: AsymmetricLossFunction ; end
type HuberLoss <: SymmetricLossFunction
	p::Float64
	HuberLoss(p::Float64) = (p <= 0) ? error("Huber loss parameter must be strictly positive") : new(p)
end
HuberLoss() = HuberLoss(1.0)
HuberLoss(p::Number) = HuberLoss(convert(Float64, p))

#----------------------------------------------------------
#LOSS FUNCTION ADMIN METHODS
#----------------------------------------------------------
#string
string(lF::SquaredLoss) = "squaredLoss"
string(lF::AbsoluteLoss) = "absoluteLoss"
string(lF::MinkowskiLoss) = "minkowskiLoss"
string(lF::SquaredLogLoss) = "squaredLogLoss"
string(lF::AbsoluteLogLoss) = "absoluteLogLoss"
string(lF::SquaredPropLoss) = "squaredProportionalLoss"
string(lF::AbsolutePropLoss) = "absoluteProportionalLoss"
string(lF::QLIKE) = "QLIKE"
string(lF::QLIKEReverse) = "QLIKEReverse"
string(lF::HuberLoss) = "huberLoss"
#show
function show{T<:Union{MinkowskiLoss, HuberLoss}}(io::IO, lf::T)
	println(io, "Loss function = " * string(lF))
	println(io, "    p = " * string(lF.p))
end
show(io::IO, lf::LossFunction) = println(io, "Loss function = " * string(lF))
show(lf::LossFunction) = show(STDOUT, lf)
#tolossfunction    WARNING: NOT TYPE STABLE
function tolossfunction(s::Symbol)
	if s == :squaredLoss ; return(SquaredLoss())
	elseif s == :absoluteLoss ; return(AbsoluteLoss())
	elseif s == :minkowskiLoss ; return(MinkowskiLoss())
	elseif s == :squaredLogLoss ; return(SquaredLogLoss())
	elseif s == :absoluteLogLoss ; return(AbsoluteLogLoss())
	elseif s == :squaredPropLoss ; return(SquaredPropLoss())
	elseif s == :absolutePropLoss ; return(AbsolutePropLoss())
	elseif s == :QLIKE ; return(QLIKE())
	elseif s == :QLIKEReverse ; return(QLIKEReverse())
	elseif s == :huberLoss ; return(HuberLoss())
	else ; error("Invalid input symbol")
	end
end
tolossfunction(x::String) = tolossfunction(symbol(x))

#----------------------------------------------------------
# CORE METHODS FOR COMPUTING LOSS
#----------------------------------------------------------
#Methods that can be defined in terms of e = x - y, i.e. single input loss methods
loss(e::Number, lF::SquaredLoss) = e^2
loss(e::Number, lF::AbsoluteLoss) = abs(e)
loss(e::Number, lF::MinkowskiLoss) = abs(e)^lF.p
loss(e::Number, lF::HuberLoss) = abs(e) <= lF.p ? 0.5*e^2 : lF.p*(abs(e)-0.5*lF.p)
#Methods that cannot be defined in terms of e = x - y
loss{T<:singleInputLossMethods}(x::Number, y::Number, lF::T) = loss(x - y)
loss(e::Number, lF::SquaredLogLoss) = (log(x) - log(y))^2
loss(e::Number, lF::AbsoluteLogLoss) = abs(log(x) - log(y))
loss(e::Number, lF::SquaredPropLoss) = (x/y - 1)^2
loss(e::Number, lF::AbsolutePropLoss) = abs(x/y - 1)
loss(e::Number, lF::QLIKE) = x/y - log(x/y) - 1
loss(e::Number, lF::QLIKEReverse) = y/x - log(y/x) - 1
#Method for a loss differential (used a lot in forecast evaluation), i.e. L(x, basecase) - L(y, basecase)
lossdiff{T<:Number}(x1::T, x2::T, xBase::T, lF::LossFunction) = loss(x1, xBase, lF) - loss(x2, xBase, lF)

#----------------------------------------------------------
# ARRAY BASED EXTENSIONS TO CORE METHODS
#----------------------------------------------------------
loss{T1<:Number, T2<:singleInputLossMethods}(e::Vector{T}, lF::LossFunction) = T[ loss(e[n], lF) for n = 1:length(e) ]
loss{T<:Number}(x::Vector{T}, y::Vector{T}, lF::LossFunction) = length(x) == length(y) ? T[ loss(x[n], y[n], lF) for n = 1:length(x) ] : error("Length mismatch")
loss{T<:Number}(x::Matrix{T}, y::Matrix{T}, lF::LossFunction) = size(x) == size(y) ? T[ loss(x[n, m], y[n, m], lF) for n = 1:size(x, 1), m = 1:size(x, 2) ] : error("size mismatch")
loss{T<:Number}(x::Vector{T}, y::Matrix{T}, lF::LossFunction) = length(x) == size(y, 1) ? T[ loss(x[n], y[n, m], lF) for n = 1:size(y, 1), m = 1:size(y, 2) ] : error("size mismatch")
loss{T<:Number}(x::Matrix{T}, y::Vector{T}, lF::LossFunction) = size(x, 1) == length(y) ? T[ loss(x[n, m], y[n], lF) for n = 1:size(x, 1), m = 1:size(x, 2) ] : error("Size mismatch")
lossdiff{T<:Number}(x1::Vector{T}, x2::Vector{T}, xBase::Vector{T}, lF::LossFunction) = (length(xBase) == length(x1) == length(x2) ? T[ lossdiff(x1[n], x2[n], xBase[n], lF) for n = 1:length(xBase) ] : error("Length mismatch")
lossdiff{T<:Number}(x1::Matrix{T}, x2::Vector{T}, xBase::Vector{T}, lF::LossFunction) = (length(xBase) == size(x1, 1) == length(x2)) ? T[ lossdiff(x1[n, m], x2[n], xBase[n], lF) for n = 1:size(x1, 1), m = 1:size(x1, 2) ]
lossdiff{T<:Number}(x1::Vector{T}, x2::Matrix{T}, xBase::Vector{T}, lF::LossFunction) = -1 * lossdiff(x2, x1, xBase, lF)

end # module


text:
using LOLTools.LoLStatusV3

api_key = get(ENV, "RIOT_TOKEN", "")
platform = "kr"

status = LoLStatusV3.shard_data(api_key, platform)
@info status


text:
#
# front-end for StructJuMP solvers interfaces
#

module StructJuMPSolverInterface

import MPI

# Struct Model interface
abstract type ModelInterface end
#abstract ModelInterface

export ModelInterface, KnownSolvers, sj_solve, getModel, getVarValue, getVarValues, getNumVars,
        getNumCons, getTotalNumVars, getTotalNumCons, getLocalBlocksIds, getLocalChildrenIds,
        getObjectiveVal

const KnownSolvers = Dict{AbstractString,Function}();

# In general we will use same ret code symbol for Ipopt
# also with added return code for PIPS
# From Ipopt/src/Interfaces/IpReturnCodes_inc.h
ApplicationReturnStatus = Dict(
  0=>:Solve_Succeeded,
  1=>:Solved_To_Acceptable_Level,
  2=>:Infeasible_Problem_Detected,
  3=>:Search_Direction_Becomes_Too_Small,
  4=>:Diverging_Iterates,
  5=>:User_Requested_Stop,
  6=>:Feasible_Point_Found,
  #for PIPS specific
  7=>:Need_Feasibility_Restoration,
  8=>:Unknown,
  #end for PIPS retcode
  -1=>:Maximum_Iterations_Exceeded,
  -2=>:Restoration_Failed,
  -3=>:Error_In_Step_Computation,
  -4=>:Maximum_CpuTime_Exceeded,
  -10=>:Not_Enough_Degrees_Of_Freedom,
  -11=>:Invalid_Problem_Definition,
  -12=>:Invalid_Option,
  -13=>:Invalid_Number_Detected,
  -100=>:Unrecoverable_Exception,
  -101=>:NonIpopt_Exception_Thrown,
  -102=>:Insufficient_Memory,
  -199=>:Internal_Error
  )

function sj_solve(model; solver="Unknown", with_prof=false, suppress_warmings=false,kwargs...)
    if !haskey(KnownSolvers,solver)
        Base.warn("Unknow solver: ", solver)
        Base.error("Known solvers are: ", keys(KnownSolvers))
    end
    status = KnownSolvers[solver](model; with_prof=with_prof, suppress_warmings=false,kwargs...)

    if !haskey(ApplicationReturnStatus,status)
      Base.warn("solver can't solve the problem");
      return :Error
    else
      # @show ApplicationReturnStatus[status]
      return ApplicationReturnStatus[status]
    end
end

# package code goes here
include("helper.jl")
# include("structure_helper.jl")
# include("nonstruct_helper.jl")

function getModel(m,id)
    return id==0?m:getchildren(m)[id]
end

function getVarValues(m,id)
    mm = getModel(m,id);
    v = Float64[];
    for i = 1:getNumVars(m,id)
        v = [v;JuMP.getvalue(JuMP.Variable(mm,i))]
    end
    return v
end

function getVarValue(m,id,idx)
    mm = getModel(m,id)
    @assert idx<=getNumVars(m,id)
    return JuMP.getvalue(JuMP.Variable(mm,idx))
end

##!feng function getObjectiveValue(m)
##! It seems that getobjectivevalue(m) does not return the correct objetive value
##!feng end
function getObjectiveVal(m)
  #MPI is already finalized , therefore returns only objective value at local nodes
  # mid, nprocs = getMyRank()
  # lobj = 0.0
  # x0 = getVarValues(m,0)

  # if mid == 0
  #   e =  get_nlp_evaluator(m,0)
  #   lobj = MathProgBase.eval_f(e,build_x(m,0,x0,x0))
  # end

  # for i in getLocalChildrenIds(m)
  #     x1 = getVarValues(m,i)
  #     lobj += MathProgBase.eval_f(get_nlp_evaluator(m,i),build_x(m,i,x0,x1))
  # end

  # obj = MPI.Reduce(lobj,MPI.SUM,0,getStructure(m).comm)
  lobj =0.0
  x0 = getVarValues(m,0)
  # @show x0
  e = get_nlp_evaluator(m,0)
  lobj = MathProgBase.eval_f(e,build_x(m,0,x0,x0))
  for i in getchildren(m)
    id = i[1]
    e = get_nlp_evaluator(m,id)
    x1 = getVarValues(m,id)
    lobj += MathProgBase.eval_f(e,build_x(m,id,x0,x1))
  end
  return lobj;
end

function getNumVars(m,id)
  mm = getModel(m,id)
  nvar = MathProgBase.numvar(mm) - length(getStructure(mm).othermap)
  return nvar
end

function getNumCons(m,id)
  mm = getModel(m,id)
  return MathProgBase.numconstr(mm)
end

function getTotalNumVars(m)
    nvar = 0
    for i=0:num_scenarios(m)
        nvar += getNumVars(m,i)
    end
    return nvar
end

function getTotalNumCons(m)
    ncon = 0
    for i=0:num_scenarios(m)
        ncon += getNumCons(m,i)
    end
    return ncon
end

function getLocalBlocksIds(m)
  myrank,mysize = getMyRank()
  numScens = num_scenarios(m)
  d = div(numScens,mysize)
  s = myrank * d + 1
  e = myrank == (mysize-1)? numScens:s+d-1
  ids=[0;s:e]
end

function getLocalChildrenIds(m)
    myrank,mysize = getMyRank()
    # @show myrank, mysize
    numScens = num_scenarios(m)
    d = div(numScens,mysize)
    s = myrank * d + 1
    e = myrank == (mysize-1)? numScens:s+d-1
    ids = collect(s:e)
end

end # module StructJuMPSolverInterface

include("pips_parallel.jl")
include("pips_serial.jl")
include("ipopt_serial.jl")


text:
mutable struct pose_node_t <: LCMType
  utime::Int64   # unix time associated with this pose
  id::Int64      # pose id (sequential, 0-based)
  mean_dim::Int32       # default will be 7 (3+4)
  mean::Vector{Float64} # default will be [x y z qw qx qy qz]
  covar_dim::Int32      # usually x y z yaw pitch roll (RotXYZ)
  covar::Vector{Float64}
end

@lcmtypesetup(pose_node_t,
  mean => (mean_dim,),
  covar => (covar_dim,)
)


text:
export Model2Fit_Empty

"""

An empty model without parameters. This is useful to initialize
summation:

```julia
model = Model2Fit_Empty()

for i in 1:n
    model = model + GaussianPeak()
end
```
"""
struct Model2Fit_Empty <: Abstract_Model2Fit
end

parameter_size(m::Model2Fit_Empty) = 0
accumulate_y!(m::Model2Fit_Empty,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = Y


text:
"""
Shunt constraint using linearized voltage magnitude difference phi

```
sum(p) + sum(pconv_grid_ac)  == sum(pg) - sum(pd) - sum(gs*(1.0 + 2*phi)
```
"""
function constraint_kcl_shunt(pm::AbstractLPACModel, n::Int, cnd::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs)
    phi = PowerModels.var(pm, n, cnd, :phi, i)
    p = PowerModels.var(pm, n, cnd, :p)
    q = PowerModels.var(pm, n, cnd, :q)
    pg = PowerModels.var(pm, n, cnd, :pg)
    qg = PowerModels.var(pm, n, cnd, :qg)
    pconv_grid_ac = PowerModels.var(pm, n, cnd, :pconv_tf_fr)
    qconv_grid_ac = PowerModels.var(pm, n, cnd, :qconv_tf_fr)

    PowerModels.con(pm, n, cnd, :kcl_p)[i] = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac)  == sum(pg[g] for g in bus_gens)   - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*(1.0 + 2*phi))
    PowerModels.con(pm, n, cnd, :kcl_q)[i] = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac)  == sum(qg[g] for g in bus_gens)   - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*(1.0 + 2*phi))
end
"""
DC branch power flow using linearized voltage magnitude difference phi

```
p_dc_fr == p * g *  (phi_fr - phi_to)
p_dc_to == p * g *  (phi_to - phi_fr)
```
"""
function constraint_ohms_dc_branch(pm::AbstractLPACModel, n::Int, cnd::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = PowerModels.var(pm, n, cnd, :p_dcgrid, f_idx)
    p_dc_to = PowerModels.var(pm, n, cnd, :p_dcgrid, t_idx)
    phi_fr = PowerModels.var(pm, n, cnd, :phi_vdcm, f_bus)
    phi_to = PowerModels.var(pm, n, cnd, :phi_vdcm, t_bus)
    phi_fr_ub = JuMP.UpperBoundRef(phi_to)
    phi_fr_lb = JuMP.LowerBoundRef(phi_to)

    if r == 0
        @constraint(pm.model, p_dc_fr + p_dc_to == 0)
    else
        g = 1 / r
        @constraint(pm.model, p_dc_fr == p * g *  (phi_fr - phi_to))
        @constraint(pm.model, p_dc_to == p * g *  (phi_to - phi_fr))
    end
end

function constraint_dc_branch_current(pm::AbstractLPACModel, n::Int, cnd::Int, f_bus, f_idx, ccm_max, p)

end

function add_dc_bus_voltage_setpoint(sol, pm::AbstractLPACModel)
    PowerModels.add_setpoint!(sol, pm, "busdc", "vm", :phi_vdcm, status_name="Vdc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using YLEFlow
using Test

@testset "YLEFlow.jl" begin

    # Write your tests here.
    g()
end


text:
Random.seed!(42)

@testset "Fourier shifting methods" begin

    @testset "Empty shifts" begin
        x = randn(ComplexF32, (11, 12, 13))
        @test FourierTools.shift(x, []) == x

        x = randn(Float32, (11, 12, 13))
        @test FourierTools.shift(x, []) == x
    end

    @testset "Integer shifts for complex and real arrays" begin
        x = randn(ComplexF32, (11, 12, 13))

        s = (2,2,2)
        @test FourierTools.shift(x, s) ≈ circshift(x, s)
        s = (3,2,1)
        @test FourierTools.shift(x, s) ≈ circshift(x, s)

        @test FourierTools.shift(x, (0,0,0)) == x
        x = randn(Float32, (11, 12, 13))

        s = (2,2,2)
        @test FourierTools.shift!(copy(x), s) ≈ circshift(x, s)
        s = (3,2,1)
        @test FourierTools.shift!(copy(x), s) ≈ circshift(x, s)

        @test sum(x) ≈ sum(FourierTools.shift!(copy(x), s))

    end

    @testset "Half integer shifts" begin

        x = [0.0, 1.0, 0.0, 1.0]
        xc = ComplexF32.(x)

        s = [0.5]
        @test FourierTools.shift!(copy(x), s) ≈ real(FourierTools.shift!(copy(xc), s))
        @test FourierTools.shift!(copy(x), s) ≈ real(FourierTools.shift!(copy(xc), 0.5))
        @test sum(x) ≈ sum(FourierTools.shift!(copy(x), s))

        @test sum(xc) ≈ sum(FourierTools.shift!(copy(xc), s))
    end

    @testset "Check shifts with soft_fraction" begin
        a = shift(delta((255,255)), (1.5,1.25),soft_fraction=0.1);
        @test abs(sum(a[real(a).<0])) < 3.0
        a = shift(delta((255,255)), (1.5,1.25),soft_fraction=0.0);
        @test abs(sum(a[real(a).<0])) > 5.0
    end

    @testset "Random shifts consistency between both methods" begin
        x = randn((11, 12, 13))
        s = randn((3,)) .* 10
        @test sum(x) ≈ sum(FourierTools.shift!(copy(x), s))
        @test FourierTools.shift!(copy(x), s) ≈ real(FourierTools.shift!(copy(x) .+ 0im, s))
        x = randn((11, 12, 13))
        s = randn((3,)) .* 10
        @test FourierTools.shift!(copy(x), s) ≈ real(FourierTools.shift!(copy(x) .+ 0im, s))
        @test sum(x) ≈ sum(FourierTools.shift!(copy(x), s))
    end

    @testset "Check revertibility for complex and real data" begin
        @testset "Complex data" begin
            x = randn(ComplexF32, (11, 12, 13))
            s = (-1.1, 12.123, 0.21)
            @test x ≈ shift(shift(x, s), .- s, fix_nyquist_frequency=true)
        end
        @testset "Real data" begin
            x = randn(Float32, (11, 12, 13))
            s = (-1.1, 12.123, 0.21)
            @test x ≈ shift(shift(x, s), .- s, fix_nyquist_frequency=true)
        end
    end

end


text:
# ---
# title: Non-local Mean Filter
# id: nlmf-demo
# cover: assets/nlmean_cover.png
# date: 2019-12-22
# author: Johnny Chen
# description: Use Non-local Mean filter to reduce gaussian noise
# ---

using ImageNoise
using TestImages, ImageShow, ImageCore, ImageQualityIndexes, ImageTransformations
using FileIO, Random #src

# First, load an image and and add gaussian noise to it

gray_img = float.(imresize(testimage("cameraman"), ratio=0.5))
n = AdditiveWhiteGaussianNoise(0.1)
noisy_img = apply_noise(gray_img, n)

# Then calling the standard `reduce_noise` API

f_nlmean = NonlocalMean(0.1)
denoised_img = reduce_noise(noisy_img, f_nlmean)

mosaicview(gray_img, noisy_img, denoised_img; nrow=1)

# Get the PSNR using ImageQualityIndexes package:

assess_psnr(gray_img, denoised_img)

mkpath("assets") #src
save("assets/nlmean_cover.png", denoised_img) #src


text:
function enablemissing(val = missing)
    Base.eval(Distributions, :(pdf(::Distribution, ::Missing) = $val))
    Base.eval(Distributions, :(logpdf(::Distribution, ::Missing) = log($val)))
    Base.eval(Distributions, :(zval(::Normal, ::Missing) = log($val)))
    nothing
end

function disablemissing()
    try
        Base.delete_method(@which pdf(Normal(), missing))
        Base.delete_method(@which logpdf(Normal(), missing))
        Base.delete_method(@which zval(Normal(), missing))
    catch
    end
end

max_stars_repo_name: maxmouchet/Statistikk

max_issues_repo_name: maxmouchet/Statistikk

max_forks_repo_name: maxmouchet/Statistikk