text:
using DelimitedFiles
using PyPlot
using RAFF

"""

    draw_problem(M; raff_output=nothing, model_str="logistic",
                 datafile="/tmp/output.txt")

    draw_problem(datafile::String="/tmp/output.txt"; kwargs...)

Draw the problem data given by a (`m`x`3`) `M` matrix. By default it
is assumed that the model is given by the logistic model. If no
arguments are given it is assumed that the data file is given in file
`/tmp/output.txt`.

If a [RAFFOutput](@ref) object is provided, then it plots the model
found and also the outliers (true and false positives).

Optional arguments:

  - `raff_output`: [RAFFOutput](@ref) object with the solution
    obtained
  - `model_str`: a string with the name of the model to be used to
    plot the solution. See [model_list](@ref) for details.

"""
function draw_problem(M; raff_output=nothing, model_str="logistic",
                      θsol=nothing)

    x = M[:, 1]
    y = M[:, 2]
    co = M[:, 3]

    true_outliers = findall(co .!= 0.0)

    PyPlot.scatter(x[co .== 0.0], y[co .== 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="o", s=50.0, linewidths=0.2)

    PyPlot.scatter(x[co .!= 0.0], y[co .!= 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="^", s=50.0, linewidths=0.2, label="Outliers")

    if raff_output != nothing

        n, model, modelstr = RAFF.model_list[model_str]

        modl1 = (x) -> model(x, raff_output.solution)

        t = minimum(x):0.01:maximum(x)
        PyPlot.plot(t, modl1.(t), color=PyPlot.cm."Set1"(2.0/9.0))

        # Draw outliers found by RAFF

        true_positives = intersect(true_outliers, raff_output.outliers)
        false_positives = setdiff(raff_output.outliers, true_positives)

        PyPlot.scatter(x[false_positives], y[false_positives],
                       color=PyPlot.cm."Pastel1"(0.0/9.0), marker="o",
                       linewidths=0.2, edgecolors="k", s=50.0, label="False positives")

        PyPlot.scatter(x[true_positives], y[true_positives],
                       color=PyPlot.cm."Pastel1"(0.0/9.0), marker="^",
                       s=50.0, linewidths=0.2, edgecolors="k", label="Identified outliers")

    end

    # Plot the true solution, if available
    if θsol != nothing

        modl2 = (x) -> model(x, θsol)

        PyPlot.plot(t, modl2.(t), color=PyPlot.cm."Set1"(1.0/9.0),
                    linestyle="--")

    end

    PyPlot.legend(loc="best")

    PyPlot.show()

    PyPlot.savefig("/tmp/figure.png", DPI=150)

end

function draw_problem(datafile::String="/tmp/output.txt"; kwargs...)

    fp = open(datafile, "r")

    N = parse(Int, readline(fp))

    M = readdlm(fp)

    close(fp)

    draw_problem(M; kwargs...)

end

function draw_comparison(M::Array{Float64, 2}, model_str, algorithms...)

    x = M[:, 1]
    y = M[:, 2]
    co = M[:, 3]

    true_outliers = findall(co .!= 0.0)

    PyPlot.scatter(x[co .== 0.0], y[co .== 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="o", s=50.0, linewidths=0.2)

    PyPlot.scatter(x[co .!= 0.0], y[co .!= 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="^", s=50.0, linewidths=0.2, label="Outliers")

    t = minimum(x):0.01:maximum(x)

    n, model, modelstr = RAFF.model_list[model_str]

    for (i, alg) in enumerate(algorithms)

        name, line, sol = alg

        modl1 = (x) -> model(x, sol)

        PyPlot.plot(t, modl1.(t), color=PyPlot.cm."Set1"(i/9.0),
                    linestyle=line, label=name)

        PyPlot.legend(loc="best")

        PyPlot.show()

        # PyPlot.savefig(outimage, DPI=600, bbox_inches="tight")

    end

end

function draw_comparison(model_str::String, algorithms...)

    fp = open("/tmp/output.txt", "r")

    N = parse(Int, readline(fp))

    M = readdlm(fp)

    close(fp)

    draw_comparison(M, model_str, algorithms...)

end

"""

Draw the points generated by the previous function.

"""
function draw_circle(data, outliers)

    np, = size(data)

    c = zeros(np)

    c[outliers] .= 1.0

    PyPlot.scatter(data[:, 1], data[:, 2], c=c, marker="o", s=50.0, linewidths=0.2,
                   cmap=PyPlot.cm["Paired"], alpha=0.9)

    PyPlot.axis("scaled")

    PyPlot.xticks([])

    PyPlot.yticks([])

    PyPlot.savefig("/tmp/circle.png", dpi=150, bbox_inches="tight")

end

"""

Draw the points and the solutions obtained. Save the picture in a file.

"""
function draw_circle_sol(M; model_str="circle", raff_output=nothing, other_sols...)

    x  = M[:, 1]
    y  = M[:, 2]
    co = M[:, 4]

    t = 0:0.1:2.1 * π

    ptx = (α, ρ, d) -> ρ * cos(α) + d[1]
    pty = (α, ρ, d) -> ρ * sin(α) + d[2]

    # Plot data

    true_outliers = findall(co .!= 0.0)

    PyPlot.scatter(x[co .== 0.0], y[co .== 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="o", s=50.0, linewidths=0.2)

    PyPlot.scatter(x[co .!= 0.0], y[co .!= 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="^", s=25.0, linewidths=0.2, label="Outliers")

    if raff_output != nothing

        n, model, modelstr = RAFF.model_list[model_str]

        fSol = raff_output.solution

        modl1x = (α) -> ptx(α, fSol[3], fSol[1:2])
        modl1y = (α) -> pty(α, fSol[3], fSol[1:2])

        PyPlot.plot(modl1x.(t), modl1y.(t), color=PyPlot.cm."Set1"(2.0/9.0))

        # Draw outliers found by RAFF

        true_positives = intersect(true_outliers, raff_output.outliers)
        false_positives = setdiff(raff_output.outliers, true_positives)

        if length(false_positives) > 0

            PyPlot.scatter(x[false_positives], y[false_positives],
                           color=PyPlot.cm."Pastel1"(0.0/9.0), marker="o",
                           linewidths=0.2, edgecolors="k", s=50.0, label="False positives")

        end

        if length(true_positives) > 0

            PyPlot.scatter(x[true_positives], y[true_positives],
                           color=PyPlot.cm."Pastel1"(0.0/9.0), marker="^",
                           s=50.0, linewidths=0.2, edgecolors="k", label="Identified outliers")

        end

    end

    PyPlot.legend(loc="best")

    PyPlot.axis("scaled")

    PyPlot.xticks([])

    PyPlot.yticks([])

    PyPlot.savefig("/tmp/circle.png", dpi=150, bbox_inches="tight")

end

function draw_circle_sol(tSol, fSol, lsSol)

    datafile = "/tmp/output.txt"

    fp = open(datafile, "r")

    N = parse(Int, readline(fp))

    M = readdlm(fp)

    close(fp)

    x  = M[:, 1]
    y  = M[:, 2]
    ρ  = M[:, 3]
    co = M[:, 4]

    t = [0:0.1:2.1 * π;]

    ptx = (α, ρ, d) -> ρ * cos(α) + d[1]
    pty = (α, ρ, d) -> ρ * sin(α) + d[2]

    # True solution

    pptx = (α) -> ptx(α, tSol[3], tSol[1:2])
    ppty = (α) -> pty(α, tSol[3], tSol[1:2])

    PyPlot.plot(pptx.(t), ppty.(t), "b--", label="True solution")

    # RAFF solution

    pptx = (α) -> ptx(α, fSol[3], fSol[1:2])
    ppty = (α) -> pty(α, fSol[3], fSol[1:2])

    PyPlot.plot(pptx.(t), ppty.(t), "g-", label="RAFF")

    # # LS solution

    # pptx = (α) -> ptx(α, lsSol[3], lsSol[1:2])
    # ppty = (α) -> pty(α, lsSol[3], lsSol[1:2])

    # PyPlot.plot(pptx.(t), ppty.(t), "r-", label="Least squares")

    PyPlot.scatter(x[co .== 0.0], y[co .== 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="o", s=50.0, linewidths=0.2)

    PyPlot.scatter(x[co .!= 0.0], y[co .!= 0.0], color=PyPlot.cm."Pastel1"(1.0/9.0),
                   marker=".", s=25.0, linewidths=0.2, label="Outliers")

    PyPlot.legend(loc=4)

    PyPlot.axis("scaled")

    PyPlot.xticks([])

    PyPlot.yticks([])

    PyPlot.savefig("/tmp/circle.png", dpi=150, bbox_inches="tight")

end

"""

    draw_img_sol(imgfile, sol; thck::Int=2)

Draw the image and the solution found, using JuliaImage
packages. `imagefile` is the path to the image and `sol` is a
3-dimensional vector with the solution.

Extra parameter `thck` defines the thickness of the circle.

"""
function draw_img_sol(imgfile, sol; thck::Int=2)

    img = load(imgfile)

    h, w = size(img)

    t = [0:0.01:2.1 * π;]

    ptx = (α, ρ, d) -> ρ * cos(α) + d[1]
    pty = (α, ρ, d) -> ρ * sin(α) + d[2]

    # Solution

    pptx = (α) -> ptx(α, sol[3], sol[1:2])
    ppty = (α) -> pty(α, sol[3], sol[1:2])

    for α in t

        i = Int(round(pptx(α)))
        j = Int(round(ppty(α)))

        !((1 <= i <= h) && (1 <= j <= w)) && continue

        img[max(1, i - thck):min(i + thck, h), max(1, j - thck):min(j + thck, w)] .= RGB(1.0, 0.0, 0.0)

    end

    ImageView.imshow(img)

end


text:
#=
Created on June 22, 2020
@author: Zhou Xiang
=#

using Images

function pad_to_same_size(complex_amp::Array{Complex{Float64}, 2},
                          lower::Int64,
                          upper::Int64)
    new_complex_amp = zeros(ComplexF64, lower + upper, lower + upper)
    new_complex_amp[lower + 1:upper, lower + 1:upper] = complex_amp
    return new_complex_amp
end

function resize_to_same_size(complex_amp_part::Array{Float64,2},
                             median_N::Int64)
    new_size = (median_N, median_N)
    resized_complex_amp_part = imresize(complex_amp_part, new_size)
    return resized_complex_amp_part
end

ext: jl

max_stars_repo_name: Miyoshichi/SimLight

max_issues_repo_name: Miyoshichi/SimLight

max_forks_repo_name: Miyoshichi/SimLight


text:
using YaoExtensions, Yao
using Test, Random
using Optim: LBFGS, optimize
using Optim

# port the `Matrix` function to Yao's AD.
using Zygote
include("chainrules_patch.jl")

function loss(u, ansatz)
    m = Matrix(ansatz)
    sum(abs.(u .- m))
end

"""
    learn_u4(u::AbstractMatrix; niter=100)

Learn a general U4 gate. The optimizer is LBFGS.
"""
function learn_u4(u::AbstractMatrix; niter=100)
    ansatz = general_U4() * put(2, 1=>phase(0.0))  # initial values are 0, here, we attach a global phase.
    params = parameters(ansatz)
    g!(G, x) = (ansatz=dispatch(ansatz, x); G .= Zygote.gradient(ansatz->loss(u, ansatz), ansatz)[1])
    optimize(x->(ansatz=dispatch(ansatz, x); loss(u, ansatz)), g!, parameters(ansatz),
                    LBFGS(), Optim.Options(iterations=niter))
    println("final loss = $(loss(u,ansatz))")
    return ansatz
end

using Random
Random.seed!(3)
u = rand_unitary(4)
c = learn_u4(u; niter=150)


text:
using DiffEqFlux, Flux
abstract type AffineSystem <: Ode end

""" dxdt = Ax + Bu + ν
    AffineOde """
struct AffineOde{FU,F2} <: AffineSystem
    state_dim::Int
    input_dim::Int
    paramsum::Int
    index_dict::Dict{String,UnitRange{Int}}
    σ_u::FU
    initial_params::F2
    function AffineOde(
        state_dim::Integer,
        input_dim::Integer;
        initW = Flux.glorot_uniform,
        initb = Flux.zeros,
        σ_u = identity,
    )
        # compute indices
        i = 1
        index_dict = Dict{String,UnitRange{Int}}()
        index_dict["A"] = i:i+state_dim*state_dim-1
        i += state_dim * state_dim
        index_dict["B"] = i:i+input_dim*state_dim-1
        i += input_dim * state_dim
        index_dict["ν"] = i:i+state_dim-1
        i += state_dim

        # initialize initial parameters
        initial_params() = vcat(
            vec(initW(state_dim, state_dim)),
            vec(initW(state_dim, input_dim)),
            vec(initb(state_dim)),
        )

        new{typeof(σ_u),typeof(initial_params)}(
            state_dim,
            input_dim,
            i - 1,
            index_dict,
            σ_u,
            initial_params,
        )
    end
end

function unpack_params(f::AffineOde, p)
    A = reshape(view(p, f.index_dict["A"]), f.state_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.state_dim, f.input_dim)
    ν = view(p, f.index_dict["ν"])
    return A, B, ν
end

function (f::AffineOde)(xu, p)
    x, u = extract_inputs(f, xu)
    A, B, ν = unpack_params(f, p)
    return A * x .+ B * f.σ_u.(u) .+ ν
end

""" dxdt = Ax + Bu
    LinearOde """
struct LinearOde{FU,F2} <: AffineSystem
    state_dim::Int
    input_dim::Int
    paramsum::Int
    index_dict::Dict{String,UnitRange{Int}}
    σ_u::FU
    initial_params::F2
    function LinearOde(
        state_dim::Integer,
        input_dim::Integer;
        initW = Flux.glorot_uniform,
        σ_u = identity,
    )
        # compute indices
        i = 1
        index_dict = Dict{String,UnitRange{Int}}()
        index_dict["A"] = i:i+state_dim*state_dim-1
        i += state_dim * state_dim
        index_dict["B"] = i:i+input_dim*state_dim-1
        i += input_dim * state_dim

        # initialize initial parameters
        initial_params() =
            vcat(vec(initW(state_dim, state_dim)), vec(initW(state_dim, input_dim)))

        new{typeof(σ_u),typeof(initial_params)}(
            state_dim,
            input_dim,
            i - 1,
            index_dict,
            σ_u,
            initial_params,
        )
    end
end

function unpack_params(f::LinearOde, p)
    A = reshape(view(p, f.index_dict["A"]), f.state_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.state_dim, f.input_dim)
    ν = zeros(f.state_dim)
    return A, B, ν
end

function (f::LinearOde)(xu, p)
    x, u = extract_inputs(f, xu)
    A = reshape(view(p, f.index_dict["A"]), f.state_dim, f.state_dim)
    B = reshape(view(p, f.index_dict["B"]), f.state_dim, f.input_dim)
    return A * x .+ B * f.σ_u.(u)
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
function ldlfor!(a::Matrix{Float64}, b::Vector{Float64})
  n = size(a, 1)
  for i in 1:n
    total = b[i]
    if i > 1
      for j in 1:(i-1)
        total -= a[j, i] * b[j]
      end
    end
    b[i] = total / a[i, i]
  end
end


text:
function test_meta(detector_type::Type{<:OD.Detector})
    @testset "$detector_type @meta" begin
        @testset "the detector `metadata_pkg` is defined" begin
            @test MLJ.package_name(detector_type) != "unknown"
            @test MLJ.package_uuid(detector_type) != "unknown"
            @test MLJ.package_url(detector_type) != "unknown"
        end

        @testset "the detector `load_path` is defined on the module top-level" begin
            @test MLJ.load_path(detector_type) == "$detector_type"
        end
    end
end


text:
using HDF5
using DocOpt
using HPAT

function generate_files(path,N)
    D = 10
    p = 3
    A = rand(D+p,N)
    hdf5_file = path*"linear_regression.hdf5"
    csv_file = path*"linear_regression.csv"
    if isfile(hdf5_file)
        rm(hdf5_file)
    end
    if isfile(csv_file)
        rm(csv_file)
    end
    h5write(hdf5_file,"/points", A[p+1:end,:])
    h5write(hdf5_file,"/labels", A[1:p,:])
    writecsv(csv_file,A')
end

function main()
    doc = """generate_linear_regression.jl

generate random input for linear regression example.

Usage:
  generate_linear_regression.jl -h | --help
  generate_linear_regression.jl [--instances=<instances>] [--path=<path>]

Options:
  -h --help                  Show this screen.
  --instances=<instances>    Specify number of instances; defaults to 2000000.
  --path=<path>              Specify output path for generated files; defaults to HPAT's default data path.
"""

    arguments = docopt(doc)

    if (arguments["--instances"] != nothing)
        instances = parse(Int, arguments["--instances"])
    else
        instances = 2*10^6
    end

    if (arguments["--path"] != nothing)
        path = arguments["--path"]
    else
        path = HPAT.getDefaultDataPath()
    end

    generate_files(path, instances)
end

main()


text:
import CUDA.CUDNN: cudnnTensorDescriptor, cudnnFilterDescriptor
using AutoGrad: Value

cudnnTensorDescriptor(x::Value) = cudnnTensorDescriptor(value(x))
cudnnFilterDescriptor(x::Value) = cudnnFilterDescriptor(value(x))


text:
#Small utility library for the definition of the State Vector
import Base.size, Base.getindex, Base.setindex!, Base.IndexStyle
import Base.*, Base.+, Base.-, Base./

!isdefined(:Quaternion) ? include("Quatlib.jl") : nothing

mutable struct State <: AbstractArray{Float64,1}
    vx::Float64
    vy::Float64
    vz::Float64
    x::Float64
    y::Float64
    z::Float64
    wx::Float64
    wy::Float64
    wz::Float64
    q1::Float64
    q2::Float64
    q3::Float64
    q4::Float64
end

Base.size(::State) = (13,)
Base.IndexStyle(::Type{State}) = IndexLinear()
function Base.getindex(state::State, i::Int)
    if (i > 0 && i < 14)
        if i == 1
            return state.vx
        elseif i == 2
            return state.vy
        elseif i == 3
            return state.vz
        elseif i == 4
            return state.x
        elseif i == 5
            return state.y
        elseif i == 6
            return state.z
        elseif i == 7
            return state.wx
        elseif i == 8
            return state.wy
        elseif i == 9
            return state.wz
        elseif i == 10
            return state.q1
        elseif i == 11
            return state.q2
        elseif i == 12
            return state.q3
        elseif i == 13
            return state.q4
        end
    else
        error("index out of range")
    end
end

function setindex!(state::State, v::T, i::Int) where {T <: Real}
    if (i > 0 && i < 14)
        if i == 1
            state.vx = Float64(v)
        elseif i == 2
            state.vy = Float64(v)
        elseif i == 3
            state.vz = Float64(v)
        elseif i == 4
            state.x = Float64(v)
        elseif i == 5
            state.y = Float64(v)
        elseif i == 6
            state.z = Float64(v)
        elseif i == 7
            state.wx = Float64(v)
        elseif i == 8
            state.wy = Float64(v)
        elseif i == 9
            state.wz = Float64(v)
        elseif i == 10
            state.q1 = Float64(v)
        elseif i == 11
            state.q2 = Float64(v)
        elseif i == 12
            state.q3 = Float64(v)
        elseif i == 13
            state.q4 = Float64(v)
        end
    else
        error("index out of range")
    end
end

velocity(s::State) = s[1:3]
position(s::State) = s[4:6]
angularvelocity(s::State) = s[7:9]
orientation(s::State) = Quaternion(s[10], s[11:13])
function State(vel::Vector{T} = zeros(3),pos::Vector{S} = zeros(3),
        angvel::Vector{R} = zeros(3), quat::Quaternion = Quaternion(0., [0., 0., 0.])) where {T,S,R <: Real}
    State(vel[1], vel[2], vel[3], pos[1], pos[2], pos[3],
        angvel[1], angvel[2], angvel[3], quat[1], quat[2], quat[3], quat[4])
end


text:
import StatsBase
import StatsBase: ZScoreTransform
export ZScoreTransform

import MultivariateStats
import MultivariateStats:
    PCA, pcacov, indim, outdim, projection, principalvar, principalvars,
    tprincipalvar, tresidualvar, tvar, principalratio, transform, reconstruct

export
    pca, PCA, pcacov, indim, outdim, projection, principalvar, principalvars,
    tprincipalvar, tresidualvar, tvar, principalratio, transform, reconstruct

"""
    pca(x::WeightedCovMatrix)::Tuple{StatsBase.ZScoreTransform, MultivariateStats.PCA}

Constructs a `StatsBase.ZScoreTransform`, and a `MultivariateStats.PCA` object
from a `WeightedCovMatrix`.

# Parameters
- `cov_pca::Bool = false`: Do a covariance or correlation pca.
- `maxoutdim::Int = size(x, 1)`: Limits the maximum number of output dimensions.
- `pratio::AbstractFloat = 1.0`: Limits the number of output dimensions to a ratio of explained variance.
- `correct::Bool = false`: Correct the cov(x)/cor(x)/std(x), see the documentation of weights in `StatsBase.jl` for details.
- `weight_type::Symbol`: The type of weight for the correction, see the documentation of weights in `StatsBase.jl` for details. Ignored if `correct` is `false`.

# Example
    c = fit!(WeightedCovMatrix(), rand(100, 4), rand(100))
    t, p = pca(c)
    # note that for prediction/reconstruction with the `PCA` object, the observations are in columns!
    transform((t, p), rand(4, 100))
    reconstruct(p, t), rand(4, 100))
"""
function pca(
        x           :: WeightedCovMatrix{T};
        cov_pca     :: Bool          = false,
        maxoutdim   :: Int           = size(x, 1),
        pratio      :: AbstractFloat = 1.0,
        correct     :: Bool          = false,
        weight_type :: Symbol        = :analytic
    ) where T

    d = size(x, 1)
    c = cov_pca ?
        cov(x, corrected = correct, weight_type = weight_type) :
        cor(x, corrected = correct, weight_type = weight_type)
    m = mean(x)
    s = cov_pca ?
        T[] :
        std(x, corrected = correct, weight_type = weight_type)

    t = ZScoreTransform(d, 1, m, s)
    p = MultivariateStats.pcacov(c, T[], maxoutdim = maxoutdim, pratio = pratio)

    return t, p
end

function transform(t::Tuple{ZScoreTransform, PCA}, x)
    xz = StatsBase.transform(t[1], x)
    transform(t[2], xz)
end

function reconstruct(t::Tuple{PCA, ZScoreTransform}, x)
    xz = reconstruct(t[1], x)
    StatsBase.reconstruct(t[2], xz)
end

transform(t::Tuple{}, x) = x
transform(t::Tuple{T}, x) where T = transform(t, x)
function transform(t::Tuple, x)
    xt = transform(x[1])
    for i in 2:length(t)
        xt = transform(t[i], xt)
    end
    xt
end

reconstruct(t::Tuple{}, x) = x
reconstruct(t::Tuple{T}, x) where T = reconstruct(t, x)
function reconstruct(t::Tuple, x)
    xt = reconstruct(x[1])
    for i in 2:length(t)
        xt = reconstruct(t[i], xt)
    end
    xt
end


text:
struct LevelSetToFractions{M, T, J, C}
    mesh::M
    topology::T
    jac::J
    colors::C
end
function LevelSetToFractions(mesh::Mesh)
    topo = topology(mesh)
    sparsity = get_sparsity_pattern(topo)
    jac = float.(sparsity)
    colors = matrix_colors(jac)
    return LevelSetToFractions(mesh, topo, jac, colors)
end
function Base.show(io::IO, ls::LevelSetToFractions)
    nelems = nelements(ls.mesh)
    nnodes = nvertices(ls.mesh)
    println(io, "Level set function: node values -> element volume fractions")
    println(io, "  #nodes = $nnodes \n  #elems = $nelems")
end

function get_sparsity_pattern(topo)
    I = Int[]
    J = Int[]
    V = Bool[]
    for (ei, e) in enumerate(elements(topo))
        for vi in indices(e)
            push!(I, ei)
            push!(J, vi)
            push!(V, true)
        end
    end
    return sparse(I, J, V, nelements(topo), nvertices(topo))
end

function (ls::LevelSetToFractions)(fracs::AbstractVector, C::AbstractVector)
    return volume_fractions!(fracs, ls.topology, ls.mesh, C)
end
function (ls::LevelSetToFractions)(C::AbstractVector)
    fracs = similar(C, nelements(ls.topology))
    return volume_fractions!(fracs, ls.topology, ls.mesh, C)
end
function ChainRulesCore.rrule(ls::LevelSetToFractions, C::AbstractVector)
    forwarddiff_color_jacobian!(ls.jac, ls, C, colorvec = ls.colors)
    return ls(C), Δ -> begin
        (NoTangent(), ls.jac' * Δ)
    end
end


text:
import PyPlot
plt = PyPlot
include("./testl1.jl")
include("./plotting.jl")

p2c = InsarTimeseries.PHASE_TO_CM
p2mm = InsarTimeseries.PHASE_TO_CM * 365 * 10

load_geolist_intlist = InsarTimeseries.load_geolist_intlist
prune_cor = InsarTimeseries.prune_cor
remove_outliers = InsarTimeseries.remove_outliers
shrink_baseline = InsarTimeseries.shrink_baseline

invert(B::AbstractArray{T,2}, u::AbstractArray{T,1}) where {T<:Number} =
    InsarTimeseries.invert_pixel_l1(u, B)
# invt(B, v) = [p2mm * (B \ v) ; p2mm * invert(B, v) ]
# gives the stacked version, sbas version, and L1 version
invt(B, v) = [p2mm * (sum(v) / sum(B));
              p2mm * (B \ v);
              p2mm * invert(B, v) ]

function prune_igrams(g, i, u, ns = 3, ct = nothing, cp = nothing, shrink = false)
    g, i, u = remove_outliers(g, i, u, mean_sigma_cutoff = ns)
    # i, u = prune_cor(i, u, cor_pixel=cp, cor_thresh=ct)
    i, u = shrink ? shrink_baseline(g, i, u) : (i, u)
    return g, i, u
end
function prunesolve(
    g,
    i,
    u,
    B,
    nsigma = 3;
    cor_thresh = nothing,
    cor_pixel = nothing,
    shrink = false,
)
    constant = true
    geo, ints, unws = prune_igrams(g, i, u, nsigma, cor_thresh, cor_pixel, shrink)
    println("Num geos removed: $(length(g) - length(geo))")
    B = InsarTimeseries.prepB(geo, ints, constant)
    invt(B, unws)
end

function demo_point(
    rowcol;
    sigma = 4,
    max_temp = 800,
    show = true,
    refpoint = nothing,
    max_date = nothing,
    unwfile = "unw_stack.h5",
    input_dset = "stack_flat_shifted",
    ignore_file = "geolist_ignore.txt",
)
    geolist, intlist, igram_idxs =
        load_geolist_intlist(unwfile, ignore_file, max_temp, max_date = max_date)
    B = InsarTimeseries.build_B_matrix(geolist, intlist)
    Blin = sum(B, dims = 2)
    geolist, intlistall, igram_idxs_all =
        load_geolist_intlist(unwfile, ignore_file, 2000, max_date = max_date)
    Ball = InsarTimeseries.build_B_matrix(geolist, intlistall)
    Blinall = sum(Ball, dims = 2)

    unw_vals = get_stack_vals(unwfile, rowcol..., 1, input_dset, igram_idxs)
    # cc_vals = get_stack_vals("cc_stack.h5", rowcol..., 1, "stack", igram_idxs);

    unw_valsall = get_stack_vals(unwfile, rowcol..., 1, input_dset, igram_idxs_all)
    # cc_valsall = get_stack_vals("cc_stack.h5", rowcol..., 1, "stack", igram_idxs_all);

    if !isnothing(refpoint)
        println("re shifting to use $refpoint as reference")
        refunw = get_stack_vals(unwfile, refpoint..., 1, input_dset, igram_idxs)
        refunwall = get_stack_vals(unwfile, refpoint..., 1, input_dset, igram_idxs_all)
        unw_vals .-= refunw
        unw_valsall .-= refunwall
    end

    # println("Solving: Shrink=true, max temp = $max_temp")
    # @show prunesolve(geolist, intlist, unw_vals, Blin, sigma, shrink=true)
    println("Solving: no outlier remove, max temp = $max_temp")
    @show prunesolve(geolist, intlist, unw_vals, Blin, 1000, shrink = false)

    println("Solving: Shrink=false, max temp = $max_temp")
    @show prunesolve(geolist, intlist, unw_vals, Blin, sigma, shrink = false)

    println("Solving: Shrink=false, all igrams")
    @show prunesolve(geolist, intlistall, unw_valsall, Blinall, sigma, shrink = false)
    !show && return

    # Scatter plot of unw vals vs baseline
    plt.figure()
    plt.scatter(Blinall, p2c * unw_valsall, label = "all")
    plt.scatter(Blin, p2c * unw_vals, label = "600 days or less")
    plt.xlabel("Baseline (days")
    plt.ylabel("CM")
    plt.title("unw vs baseline")

    # plt.figure()
    # plt.scatter(Blinall, cc_valsall, label="all")
    # plt.scatter(Blin, cc_vals, label="$max_temp days or less")
    # plt.xlabel("Baseline (days")
    # plt.ylabel("CM")
    # plt.title("correlation vs baseline")

    # Plot unregularized
    plt.figure()
    unreg =
        p2c .* InsarTimeseries.integrate_velocities(
            B \ unw_vals,
            InsarTimeseries.day_diffs(geolist),
        )
    plt.plot(geolist, unreg, label = "unreg", "bo-")
    plt.title("Unregularized solution")
    plt.ylabel("CM")

    plot_grouped_by_day(geolist, intlist, unw_vals, sigma)
    plot_big_days(geolist, intlist, unw_vals, nsigma = sigma)
    return unreg
end