text:
#=

    This is a port of the R2 model Hiv.cs

    with the datafiles
    - hiv_hiv_y.csv
    - hiv_person.csv
    - hiv_time.csv
    - hiv_num_persons.csv

    The R2 model output:
    '''
    ...
    [168] Mean: 0.169296        [muA1]
    [168] Variance: 0.206532

    [169] Mean: 0.12894         [muA2]
    [169] Variance: 0.469274

    [170] Mean: 55.8414         [sigmaA1]
    [170] Variance: 307.621

    [171] Mean: 49.4251         [sigmaA2]
    [171] Variance: 32.243

    [172] Mean: 44.0546         [sigmaY]
    [172] Variance: 228.657
    '''

    Using HMC(0.1,6):
Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

        muA1   -0.4348    0.9159     0.0290    0.1427    42.9091    1.0093        0.4769
        muA2    0.2805    1.0354     0.0327    0.1446    56.6151    1.0393        0.6293
     sigmaA1   72.8204    2.6109     0.0826    0.1307   400.6981    0.9990        4.4539
     sigmaA2   29.7999    1.1906     0.0377    0.1017    88.2580    1.0322        0.9810
      sigmaY   71.2682    3.6674     0.1160    0.8280     3.7999    1.4849        0.0422

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

        muA1   -2.0894   -1.0887   -0.4788    0.2089    1.3453
        muA2   -1.6757   -0.4626    0.1953    0.9684    2.4151
     sigmaA1   67.2844   71.1223   72.8787   74.5292   77.6883
     sigmaA2   27.4074   28.9708   29.7822   30.6076   32.3003
      sigmaY   65.0158   68.4093   70.8732   74.1892   78.2249

138.288736 seconds (1.55 G allocations: 331.642 GiB, 15.40% gc time, 0.00% compilation time)

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function hiv(hiv_y,hiv_person,hiv_time, hiv_num_persons)
    n = length(hiv_y)

    muA1 ~ Normal(0,1)
    muA2 ~ Normal(0,1)

    sigmaA1 ~ Uniform(0,100)
    sigmaA2 ~ Uniform(0,100)

    a1 ~ filldist(Normal(muA1,sigmaA1),n)
    a2 ~ filldist(Normal(0.1*muA2,sigmaA2),n)

    sigmaY ~ Uniform(0,100)
    yHat = Vector{Real}(undef, n)
    for i in 1:n
        yHat[i] = a1[hiv_person[i]] + a2[hiv_person[i]] * hiv_time[i]
        hiv_y[i] ~ Normal(yHat[i], sigmaY)
    end

end

hiv_y = parse.(Float64,split(readline("hiv_y.csv"),","))
hiv_time  = parse.(Float64,split(readline("hiv_time.csv"),","))
hiv_num_persons = parse.(Int,split(readline("hiv_num_persons.csv"),","))
hiv_person = parse.(Int,split(readline("hiv_person.csv"),","))

println("len(y):$(length(hiv_y))")

model = hiv(hiv_y,hiv_person,hiv_time, hiv_num_persons[1])

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 1_000) # very slow
chns = sample(model, HMC(0.05,6), 1_000)

display(chns)

display(chns[[:muA1,:muA2,:sigmaA1,:sigmaA2,:sigmaY]])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using AtomIO
using Test

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

ext: jl


text:
# Functions for loading FITS files from the SloanDigitalSkySurvey.
module SDSSIO

import FITSIO
import WCS
using CodecZlib
using CodecBzip2

import ..Log
import ..Model: RawPSF, Image, CatalogEntry, eval_psf, SkyIntensity
import ..PSF
import Base.convert, Base.getindex
export PlainFITSStrategy, RunCamcolField

# types of things to mask in read_mask().
const DEFAULT_MASK_PLANES = ["S_MASK_INTERP",  # bad pixel (was interpolated)
                             "S_MASK_SATUR",  # saturated
                             "S_MASK_CR",  # cosmic ray
                             "S_MASK_GHOST"]  # electronics artifacts

const BAND_CHAR_TO_NUM = Dict('u'=>1, 'g'=>2, 'r'=>3, 'i'=>4, 'z'=>5)

struct RunCamcolField
    run::Int16
    camcol::UInt8
    field::Int16
end

function RunCamcolField(run::String, camcol::String, field::String)
    RunCamcolField(
        parse(Int16, run),
        parse(UInt8, camcol),
        parse(Int16, field))
end

abstract type SDSSImageDesc; end

struct PhotoObj <: SDSSImageDesc
    rcf::RunCamcolField
end
struct PhotoField <: SDSSImageDesc
    run::Int16
    camcol::UInt8
end
struct PsField <: SDSSImageDesc
    rcf::RunCamcolField
end
struct Frame <: SDSSImageDesc
    rcf::RunCamcolField
    band::Char
end
struct Mask <: SDSSImageDesc
    rcf::RunCamcolField
    band::Char
end
struct FieldExtents <: SDSSImageDesc; end
rcf(img::SDSSImageDesc) = img.rcf
rcf(img::PhotoField) = RunCamcolField(img.run, img.camcol, -1)

filename(p::PhotoObj) = "photoObj-$(dec(p.rcf.run,6))-$(dec(p.rcf.camcol))-$(dec(p.rcf.field,4)).fits"
filename(p::PhotoField) = "photoField-$(dec(p.run,6))-$(p.camcol).fits"
filename(p::PsField) = "psField-$(dec(p.rcf.run,6))-$(p.rcf.camcol)-$(dec(p.rcf.field,4)).fit"
filename(p::Frame) = "frame-$(p.band)-$(dec(p.rcf.run,6))-$(p.rcf.camcol)-$(dec(p.rcf.field,4)).fits"
filename(p::Mask) = "fpM-$(dec(p.rcf.run,6))-$(p.band)$(p.rcf.camcol)-$(dec(p.rcf.field,4)).fit"
filename(p::FieldExtents) = "field_extents.fits"

"""
read_sky(hdu)

Construct an image of the sky from the sky HDU of a SDSS \"frame\" file.
This is typically the 3rd HDU of the file.
"""
function read_sky(hdu::FITSIO.TableHDU)
    sky_small = squeeze(read(hdu, "ALLSKY"), 3)::Array{Float32, 2}
    sky_x = vec(read(hdu, "XINTERP"))::Vector{Float32}
    sky_y = vec(read(hdu, "YINTERP"))::Vector{Float32}

    # convert sky interpolation coordinates from 0-indexed to 1-indexed
    for i=1:length(sky_x)
        sky_x[i] += 1.0f0
    end
    for i=1:length(sky_y)
        sky_y[i] += 1.0f0
    end

    # sky intensity has to be strictly greater than 0 for the elbo
    # to be defined
    @assert all((x)-> x > 1e-12, sky_small)

    return sky_small, sky_x, sky_y
end

"""
read_frame(fname)

Read an SDSS \"frame\" FITS file and return a 4-tuple:

- `image`: sky-subtracted and calibrated 2-d image.
- `calibration`: 1-d \"calibration\".
- `sky`: 2-d sky that was subtracted.
- `wcs`: WCSTransform constructed from image header.
"""
function read_frame(strategy, rcf, b)
    f = readFITS(strategy, Frame(rcf, b))
    hdr = FITSIO.read_header(f[1], String)::String
    image = read(f[1])::Array{Float32, 2}  # sky-subtracted & calibrated data, in nMgy
    calibration = read(f[2])::Vector{Float32}
    sky_small, sky_x, sky_y = read_sky(f[3])
    close(f)

    sky = SkyIntensity(sky_small, sky_x, sky_y, calibration)

    wcs = WCS.from_header(hdr)[1]

    return image, calibration, sky, wcs
end

"""
read_field_gains(fname, fieldnum)

Return the image gains for field number `fieldnum` in an SDSS
\"photoField\" file `fname`.
"""
function read_field_gains(strategy, rcf)
    fieldnum = rcf.field
    f = readFITS(strategy, PhotoField(rcf.run, rcf.camcol))
    fieldnums = read(f[2], "FIELD")::Vector{Int32}
    gains = read(f[2], "GAIN")::Array{Float32, 2}
    close(f)

    # Find first occurance of `fieldnum` and return the corresponding gain.
    for i=1:length(fieldnums)
        if fieldnums[i] == fieldnum
            bandgains = gains[:, i]
            return Dict(zip("ugriz", bandgains))
        end
    end

    error("field number $fieldnum not found in file: $fname")
end

"""
read_mask(fname[, mask_planes])

Read a \"fpM\"-format SDSS file and return masked image ranges,
based on `mask_planes`. Returns two `Vector{UnitRange{Int}}`,
giving the range of x and y indicies to be masked.
"""
function read_mask(strategy, rcf, b, mask_planes=DEFAULT_MASK_PLANES; data=nothing)
    f = readFITS(strategy, Mask(rcf, b))

    # The last (12th) HDU contains a key describing what each of the
    # other HDUs are. Use this to find the indicies of all the relevant
    # HDUs (those with attributeName matching a value in `mask_planes`).
    value = read(f[12], "Value")::Vector{Int32}
    def = read(f[12], "defName")::Vector{String}
    attribute = read(f[12], "attributeName")::Vector{String}

    # initialize return values
    xranges = UnitRange{Int}[]
    yranges = UnitRange{Int}[]

    # Loop over keys and check if each is a mask plane we're interested in
    # (those with defName == "S_MASKTYPE" and "attributeName" in mask_planes).
    # If so, read from the corresponding HDU and construct ranges to mask.
    for i=1:length(value)
        if (def[i] == "S_MASKTYPE" && attribute[i] in mask_planes)

            # `value` starts from 0, but first table hdu is hdu number 2
            hdunum = value[i] + 2

            cmin = read(f[hdunum], "cmin")::Vector{Int32}
            cmax = read(f[hdunum], "cmax")::Vector{Int32}
            rmin = read(f[hdunum], "rmin")::Vector{Int32}
            rmax = read(f[hdunum], "rmax")::Vector{Int32}

            # "c" ("column") refers to mask's NAXIS1 (x axis);
            # "r" ("row") refers to mask's NAXIS2 (y axis).
            # cmin/cmax and rmin/rmax are 0-based and inclusive, so we add 1
            # to both.
            for j=1:length(cmin)
                push!(xranges, (cmin[j] + 1):(cmax[j] + 1))
                push!(yranges, (rmin[j] + 1):(rmax[j] + 1))
            end
        end
    end
    close(f)

    return xranges, yranges
end

struct RawImage
    rcf::RunCamcolField
    b::Int  # band index
    pixels::Matrix{Float32} # in nMgy, sky-subtracted and calibrated
    calibration::Vector{Float32}
    sky::SkyIntensity
    wcs::WCS.WCSTransform
    gain::Float32
    raw_psf_comp::RawPSF
end

"""
Read a SDSS run/camcol/field into a vector of RawImages
"""
function load_raw_images(strategy, rcf::RunCamcolField, drop_quickly::Bool = false)
    # read gain for each band
    gains = read_field_gains(strategy, rcf)

    # open FITS file containing PSF for each band
    psffile = readFITS(strategy, PsField(rcf))

    raw_images = Vector{RawImage}(5)

    for (b, band) in enumerate(['u', 'g', 'r', 'i', 'z'])
        # load image data
        r = parse_image(strategy, rcf, b, psffile, band, gains)
        if !drop_quickly
            raw_images[b] = r
        end
    end
    close(psffile)

    return raw_images
end

function parse_image(strategy, rcf, b, psffile, band, gains)
    # load image data
    pixels, calibration, sky, wcs = read_frame(strategy, rcf, band)

    # read mask
    mask_xranges, mask_yranges = read_mask(strategy, rcf, band)

    # apply mask
    for i=1:length(mask_xranges)
        pixels[mask_xranges[i], mask_yranges[i]] = NaN
    end

    # read the psf
    raw_psf_comp = read_psf(psffile, band)

    # Set it to use a constant background but include the non-constant data.
    r = RawImage(rcf, b,
                        pixels, calibration, sky, wcs,
                        gains[band], raw_psf_comp)
    r
end

"""
Load all the images for multiple rcfs
"""
function load_raw_images(strategy, rcfs::Vector{RunCamcolField}, drop_quickly::Bool = false)
    raw_images = RawImage[]

    for rcf in rcfs
        #Log.info("loading images for $rcf")
        rcf_raw_images = SDSSIO.load_raw_images(strategy, rcf, drop_quickly)
        append!(raw_images, rcf_raw_images)
    end

    raw_images
end

"""
Converts a raw image to an image by fitting the PSF, a mixture of Gaussians.
"""
function convert(::Type{Image}, r::RawImage)
    H, W = size(r.pixels)

    psfstamp = eval_psf(r.raw_psf_comp, H / 2., W / 2.)
    celeste_psf = PSF.fit_raw_psf_for_celeste(psfstamp, 2)[1]

    # scale pixels to raw electron counts
    @assert size(r.pixels, 1) == length(r.calibration)
    @inbounds for j=1:size(r.pixels, 2), i=1:size(r.pixels, 1)
        sky_ij = r.sky[i, j]
        r.pixels[i, j] = (r.gain / r.calibration[i]) * (r.pixels[i, j] + sky_ij)
    end

    nelec_per_nmgy = r.gain ./ r.calibration

    Image(H, W, r.pixels, r.b, r.wcs, celeste_psf,
          r.rcf.run, r.rcf.camcol, r.rcf.field,
          r.sky, nelec_per_nmgy, r.raw_psf_comp)
end

# -----------------------------------------------------------------------------
# PSF-related functions

"""
read_psf(fitsfile, band)

Read PSF components for the given band ('u', 'g', 'r', 'i', or 'z')
from the open FITSIO.FITS instance `fitsfile`, which should be a SDSS
\"psField\" file, and return a RawPSF instance representing the spatially
variable PSF.
"""
function read_psf(fitsfile::FITSIO.FITS, band::Char)

    # get the HDU, assuming PSF hdus are in order after primary header
    extnum = 1 + BAND_CHAR_TO_NUM[band]
    hdu = fitsfile[extnum]::FITSIO.TableHDU

    nrows = FITSIO.read_key(hdu, "NAXIS2")[1]::Int
    nrow_b = Int((read(hdu, "nrow_b")::Vector{Int32})[1])
    ncol_b = Int((read(hdu, "ncol_b")::Vector{Int32})[1])
    rnrow = (read(hdu, "rnrow")::Vector{Int32})[1]
    rncol = (read(hdu, "rncol")::Vector{Int32})[1]
    cmat_raw = read(hdu, "c")::Array{Float32, 3}
    rrows_raw = read(hdu, "rrows")::Array{Array{Float32,1},1}

    # Only the first (nrow_b, ncol_b) submatrix of cmat is used for
    # reasons obscure to the author.
    cmat = Array{Float64,3}(nrow_b, ncol_b, size(cmat_raw, 3))
    for k=1:size(cmat_raw, 3), j=1:nrow_b, i=1:ncol_b
        cmat[i, j, k] = cmat_raw[i, j, k]
    end

    # convert rrows to Array{Float64, 2}, assuming each row is the same length.
    rrows = Matrix{Float64}(length(rrows_raw[1]), length(rrows_raw))
    for i=1:length(rrows_raw)
        rrows[:, i] = rrows_raw[i]
    end

    return RawPSF(rrows, rnrow, rncol, cmat)
end

"""
Read a source catalog from an SDSS \"photoObj\" FITS file for a given
run/camcol/field combination, returning a Vector of columns and a Vector
of column names.

This is currently pretty specific to Celeste because it returns only columns
that we're interested in, and it does some transformations of some of the
columns.

The photoObj file format is documented here:
https://data.sdss.org/datamodel/files/BOSS_PHOTOOBJ/RERUN/RUN/CAMCOL/photoObj.html
"""
read_photoobj(strategy, rcf, band::Char='r') =
    read_photoobj(readFITS(strategy, PhotoObj(rcf)), band, true)

function read_photoobj(f::FITSIO.FITS, band::Char='r', close_file=true)
    b = BAND_CHAR_TO_NUM[band]

    # sometimes the expected table extension is only an empty generic FITS
    # header, indicating no objects. We check explicitly for this case here
    # and return an empty table
    if (length(f) < 2) || (!isa(f[2], FITSIO.TableHDU))
        catalog = Dict("objid"=>String[],
                       "ra"=>Float64[],
                       "dec"=>Float64[],
                       "thing_id"=>Int32[],
                       "mode"=>UInt8[],
                       "is_star"=>BitVector(),
                       "is_gal"=>BitVector(),
                       "frac_dev"=>Float32[],
                       "ab_exp"=>Float32[],
                       "theta_exp"=>Float32[],
                       "phi_exp"=>Float32[],
                       "ab_dev"=>Float32[],
                       "theta_dev"=>Float32[],
                       "phi_dev"=>Float32[],
                       "phi_offset"=>Float32[])
        for (b, n) in BAND_CHAR_TO_NUM
            catalog["psfflux_$b"] = Float32[]
            catalog["compflux_$b"] = Float32[]
            catalog["expflux_$b"] = Float32[]
            catalog["devflux_$b"] = Float32[]
        end
        return catalog
    end

    hdu = f[2]::FITSIO.TableHDU

    objid = read(hdu, "objid")::Vector{String}
    ra = read(hdu, "ra")::Vector{Float64}
    dec = read(hdu, "dec")::Vector{Float64}
    mode = read(hdu, "mode")::Vector{UInt8}
    thing_id = read(hdu, "thing_id")::Vector{Int32}

    # Get "bright" objects.
    # (In objc_flags, the bit pattern Int32(2) corresponds to bright objects.)
    objc_flags = read(hdu, "objc_flags")::Vector{Int32}
    objc_flags2 = read(hdu, "objc_flags2")::Vector{Int32}

    # bright, saturated, or large
    bad_flags1 = objc_flags .& UInt32(1^2 + 2^18 + 2^24) .!= 0

    # nopeak, DEBLEND_DEGENERATE, or saturated center
    bad_flags2 = objc_flags2 .& UInt32(2^14 + 2^18 + 2^11) .!= 0

    has_child = read(hdu, "nchild")::Vector{Int16} .> 0

    # 6 = star, 3 = galaxy, others can be ignored.
    objc_type = read(hdu, "objc_type")::Vector{Int32}
    is_star = objc_type .== 6
    is_gal = objc_type .== 3
    bad_type = ((!).(is_star .| is_gal))

    fracdev = vec((read(hdu, "fracdev")::Matrix{Float32})[b, :])

    # determine mask for rows
    has_dev = fracdev .> 0.
    has_exp = fracdev .< 1.
    is_comp = (has_dev .& has_exp)
    is_bad_fracdev = ((fracdev .< 0.) .| (fracdev .> 1))

    # TODO: We don't really want to exclude objects entirely just for being
    # bright: we just don't want to use for scoring (since
    # they're very saturated, presumably).
    mask = (!).(is_bad_fracdev .| bad_type .| bad_flags1 .| bad_flags2 .| has_child)

    # Read the fluxes.
    # Record the cmodelflux if the galaxy is composite, otherwise use
    # the flux for the appropriate type.
    psfflux = read(hdu, "psfflux")::Matrix{Float32}
    cmodelflux = read(hdu, "cmodelflux")::Matrix{Float32}
    devflux = read(hdu, "devflux")::Matrix{Float32}
    expflux = read(hdu, "expflux")::Matrix{Float32}

    # We actually only store the following properties for the given band
    # (see below in Dict construction) NB: the phi quantites in tractor are
    # multiplied by -1.
    phi_dev_deg = read(hdu, "phi_dev_deg")::Matrix{Float32}
    phi_exp_deg = read(hdu, "phi_exp_deg")::Matrix{Float32}
    phi_offset = read(hdu, "phi_offset")::Matrix{Float32}

    theta_dev = read(hdu, "theta_dev")::Matrix{Float32}
    theta_exp = read(hdu, "theta_exp")::Matrix{Float32}

    ab_exp = read(hdu, "ab_exp")::Matrix{Float32}
    ab_dev = read(hdu, "ab_dev")::Matrix{Float32}

    close_file && close(f)

    # construct result catalog
    catalog = Dict("objid"=>objid[mask],
                   "ra"=>ra[mask],
                   "dec"=>dec[mask],
                   "thing_id"=>thing_id[mask],
                   "mode"=>mode[mask],
                   "is_star"=>is_star[mask],
                   "is_gal"=>is_gal[mask],
                   "frac_dev"=>fracdev[mask],
                   "ab_exp"=>vec(ab_exp[b, mask]),
                   "theta_exp"=>vec(theta_exp[b, mask]),
                   "phi_exp"=>vec(phi_exp_deg[b, mask]),
                   "ab_dev"=>vec(ab_dev[b, mask]),
                   "theta_dev"=>vec(theta_dev[b, mask]),
                   "phi_dev"=>vec(phi_dev_deg[b, mask]),
                   "phi_offset"=>vec(phi_offset[b, mask]))
    for (b, n) in BAND_CHAR_TO_NUM
        catalog["psfflux_$b"] = vec(psfflux[n, mask])
        catalog["compflux_$b"] = vec(cmodelflux[n, mask])
        catalog["expflux_$b"] = vec(expflux[n, mask])
        catalog["devflux_$b"] = vec(devflux[n, mask])
    end

    # left over note about fluxes above: Use the cmodelflux if the
    # galaxy is composite, otherwise use the flux for the appropriate type.

    return catalog
end

"""
Convert from a catalog in dictionary-of-arrays, as returned by
read_photoobj to Vector{CatalogEntry}.
"""
function convert(::Type{Vector{CatalogEntry}}, catalog::Dict)
    out = CatalogEntry[]

    for i=1:length(catalog["objid"])
        worldcoords = [catalog["ra"][i], catalog["dec"][i]]
        frac_dev = catalog["frac_dev"][i]

        # Fill star and galaxy flux
        star_fluxes = zeros(5)
        gal_fluxes = zeros(5)
        for (j, band) in enumerate(['u', 'g', 'r', 'i', 'z'])
            # Make negative fluxes positive.
            psfflux = max(catalog["psfflux_$band"][i], 1e-6)
            devflux = max(catalog["devflux_$band"][i], 1e-6)
            expflux = max(catalog["expflux_$band"][i], 1e-6)

            # galaxy flux is a weighted sum of the two components.
            star_fluxes[j] = psfflux
            gal_fluxes[j] = frac_dev * devflux + (1.0 - frac_dev) * expflux
        end

        # For shape parameters, we use the dominant component (dev or exp)
        usedev = frac_dev > 0.5
        fits_ab = usedev? catalog["ab_dev"][i] : catalog["ab_exp"][i]
        fits_phi = usedev? catalog["phi_dev"][i] : catalog["phi_exp"][i]
        fits_theta = usedev? catalog["theta_dev"][i] : catalog["theta_exp"][i]

        # effective radius
        re_arcsec = max(fits_theta, 1. / 30)
        re_pixel = re_arcsec / 0.396

        fits_phi -= catalog["phi_offset"][i]

        # fits_phi is now degrees counter-clockwise from vertical
        # in pixel coordinates.

        # Celeste's phi measures radians counter-clockwise from vertical.
        celeste_phi_rad = fits_phi * (pi / 180)

        entry = CatalogEntry(worldcoords, catalog["is_star"][i], star_fluxes,
                             gal_fluxes, frac_dev, fits_ab, celeste_phi_rad,
                             re_pixel)
        push!(out, entry)
    end

    return out
end

function assemble_catalog(rawcatalogs::Vector{Dict};
                          duplicate_policy=:primary)
    # Limit each catalog to primary objects and objects where thing_id != -1
    # (thing_id == -1 indicates that the matching process failed)
    for cat in rawcatalogs
        mask = (cat["thing_id"] .!= -1)
        if duplicate_policy == :primary
            mask = (mask .& (cat["mode"] .== 0x01))
        end
        for key in keys(cat)
            cat[key] = cat[key][mask]
        end
    end

    #for i in eachindex(fts)
    #    Log.info(string("field $(fts[i]): $(length(rawcatalogs[i]["objid"])) ",
    #            "filtered entries"))
    #end

    # Merge all catalogs together (there should be no duplicate objects,
    # because for each object there should only be one "primary" occurance.)
    rawcatalog = deepcopy(rawcatalogs[1])
    for i=2:length(rawcatalogs)
        for key in keys(rawcatalog)
            append!(rawcatalog[key], rawcatalogs[i][key])
        end
    end

    # check that there are no duplicate thing_ids (see above comment)
    if length(Set(rawcatalog["thing_id"])) < length(rawcatalog["thing_id"])
        error("Found one or more duplicate primary thing_ids in photoobj " *
              "catalogs")
    end

    # convert to celeste format catalog
    catalog = convert(Vector{CatalogEntry}, rawcatalog)

    return catalog
end

"""
read_photoobj_files(fieldids, dirs) -> Vector{CatalogEntry}

Combine photoobj catalogs for the given overlapping fields, returning a single
joined catalog.

The `duplicate_policy` argument controls how catalogs are joined.
With `duplicate_policy = :primary`, only primary objects are included in the
combined catalog.
With `duplicate_policy = :first`, only the first detection is included in the
combined catalog.
"""
function read_photoobj_files(strategy, fts::Vector{RunCamcolField};
                             duplicate_policy=:primary,
                             slurp::Bool = false, drop_quickly::Bool = false)
    @assert duplicate_policy == :primary || duplicate_policy == :first
    @assert duplicate_policy == :primary || length(fts) == 1

    #Log.info("reading photoobj catalogs for $(length(fts)) fields")

    # the code below assumes there is at least one field.
    if length(fts) == 0
        return CatalogEntry[]
    end

    # Read in all photoobj catalogs.
    rawcatalogs = Vector{Dict}(length(fts))
    for i in eachindex(fts)
        ft = fts[i]
        #Log.info("field $(fts[i]): reading $fname")
        po = read_photoobj(strategy, ft, 'r')
        if !drop_quickly
            rawcatalogs[i] = po
        end
    end

    if drop_quickly
        return CatalogEntry[]
    end

    return assemble_catalog(rawcatalogs; duplicate_policy=duplicate_policy)
end

# IO Strategies, these are specification objects fed in from external
abstract type IOStrategy end

struct PlainFITSStrategy <: IOStrategy
    stagedir::String
    dirlayout::Symbol
    compressed::Bool
    # Could be a string or rcf -> stagedir function
    rcf_stagedir
    slurp::Bool
end
PlainFITSStrategy(stagedir::String, slurp::Bool = false) = PlainFITSStrategy(stagedir, :celeste, false, stagedir, slurp)

function slurp_fits(fname)
    is_linux() || is_bsd() || error("Slurping not implemented on this OS")
    # Do this with explicit POSIX calls, to make sure to avoid
    # any well-intended, but ultimately unhelpful intermediate
    # buffering.
    fd = ccall(:open, Cint, (Ptr{UInt8}, Cint, Cint), fname, Base.Filesystem.JL_O_RDONLY, 0)
    systemerror("open($fname)", fd == -1)
    stat_struct = zeros(UInt8, ccall(:jl_sizeof_stat, Int32, ()))
    # Can't use Base.stat because threading
    ret = ccall(:jl_fstat, Cint, (Cint, Ptr{UInt8}), fd, stat_struct)
    ret == 0 || throw(Base.UVError("stat",r))
    size = Base.Filesystem.StatStruct(stat_struct).size
    data = Array{UInt8}(size)
    rsize = ccall(:read, Cint, (Cint, Ptr{UInt8}, Csize_t), fd, data, size)
    systemerror("read", rsize != size)
    r = ccall(:close, Cint, (Cint,), fd)
    systemerror("close", r == -1)
    data
end

function compute_fname(strategy::PlainFITSStrategy, img::SDSSImageDesc)
    isa(img, FieldExtents) && return (joinpath(strategy.stagedir, filename(img)), nothing)
    imgrcf = rcf(img)
    basedir = isa(strategy.rcf_stagedir, String) ? strategy.rcf_stagedir : strategy.rcf_stagedir(imgrcf)
    if strategy.dirlayout == :celeste
        subdir = "$basedir/$(imgrcf.run)/$(imgrcf.camcol)"
        !isa(img, PhotoField) && (subdir = joinpath(subdir, string(imgrcf.field)))
    elseif strategy.dirlayout == :sdss
        if typeof(img) <: Union{Mask, PsField}
            # Photometric reductions
            subdir = joinpath(basedir, "boss/photo/redux/301", "$(imgrcf.run)/objcs/$(imgrcf.camcol)")
        elseif typeof(img) <: PhotoField
            subdir = joinpath(basedir, "boss/photoObj/301", "$(imgrcf.run)")
        elseif typeof(img) <: PhotoObj
            subdir = joinpath(basedir, "boss/photoObj/301", "$(imgrcf.run)/$(imgrcf.camcol)")
        elseif typeof(img) <: Frame
            subdir = joinpath(basedir, "boss/photoObj/frames/301", "$(imgrcf.run)/$(imgrcf.camcol)")
        end
    end
    fname = joinpath(subdir, filename(img))
    compression = nothing
    (strategy.compressed && isa(img, Frame)) && (fname *= ".bz2"; compression = Bzip2Decompression())
    (strategy.compressed && isa(img, Mask)) && (fname *= ".gz"; compression = GzipDecompression())
    fname, compression
end

function readFITS(strategy::PlainFITSStrategy, img::SDSSImageDesc)
    fname, compression = compute_fname(strategy, img)
    if strategy.slurp || compression != nothing
        data = slurp_fits(fname)
        (compression != nothing) && (data = transcode(compression, data))
        return FITSIO.FITS(data)
    else
        return FITSIO.FITS(fname)
    end
end

abstract type NetworkIoStrategy <: IOStrategy; end

struct MasterRPCStrategy <: NetworkIoStrategy
    remote_strategy::IOStrategy
end

function readRawData(strategy, img)
    fname, compression = compute_fname(strategy, img)
    data = slurp_fits(fname)
    data, compression
end

function readFITS(strategy::MasterRPCStrategy, img::SDSSImageDesc)
    data, compression = remotecall_fetch(readRawData, 1, strategy.remote_strategy, img)
    (compression != nothing) && (data = transcode(compression, data))
    return FITSIO.FITS(data)
end

"""
Read a SDSS run/camcol/field into an array of Images.
"""
function load_field_images(strategy::IOStrategy, rcfs, drop_quickly::Bool = false)
    raw_images = SDSSIO.load_raw_images(strategy, rcfs, drop_quickly)

    N = length(raw_images)
    images = Vector{Image}(N)

    drop_quickly && return images

    #Threads.@threads for n in 1:N
    for n in 1:N
        try
            images[n] = convert(Image, raw_images[n])
        catch exc
            Log.exception(exc)
            rethrow()
        end
    end

    return images
end

end  # module


text:
using RogerLuo
using Test

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


text:
using JuMP
using ECOS
using LinearAlgebra
using Random

Random.seed!(1234);

u0 = rand(10)
p = rand(10)
q = rand();

model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, u[1:10])
@variable(model, t)
@objective(model, Min, t)
@constraint(model, [t, (u - u0)...] in SecondOrderCone())
@constraint(model, u' * p == q)
optimize!(model)

@show objective_value(model);
@show value.(u);

e1 = [1, zeros(10)...]
dual_model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(dual_model, y1 <= 0)
@variable(dual_model, y2[1:11])
@objective(dual_model, Max, q * y1 + dot(vcat(0, u0), y2))
@constraint(dual_model, e1 - [0, p...] .* y1 - y2 .== 0)
@constraint(dual_model, y2 in SecondOrderCone())
optimize!(dual_model)

@show objective_value(dual_model);

model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, u[1:10])
@variable(model, t)
@objective(model, Min, t)
@constraint(model, [t, 0.5, (u - u0)...] in RotatedSecondOrderCone())
@constraint(model, u' * p == q)
optimize!(model)

@show value.(u);

n = 15;
m = 10;
A = randn(m, n);
b = rand(m, 1);

model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
@variable(model, t[1:n])
@variable(model, x[1:n])
@objective(model, Max, sum(t))
@constraint(model, sum(x) == 1)
@constraint(model, A * x .<= b )
# Cannot use the exponential cone directly in JuMP, hence we use MOI to specify the set.
@constraint(model, con[i = 1:n], [1, x[i], t[i]] in MOI.ExponentialCone())

optimize!(model);

@show objective_value(model);


text:
struct Environments
    I::MPS
    H::MPS
    InProgress::Matrix{ITensor}
    DiagInProgress::Matrix{ITensor}
end

function Base.copy(E::Environments)
    return Environments(copy(E.I), copy(E.H), copy(E.InProgress))
end

function getGuessInds(A::fPEPS, prev_mps, col::Int)
    # guess should have the indices of A going left/right that prev_mps does not have
    Ny, Nx = size(A)
    A_prev_common  = [commonind(A[row, col], prev_mps[row]) for row in 1:Ny]
    hori_A_inds    = [inds(A[row, col], "Link, r") for row in 1:Ny]
    hori_prev_inds = [inds(prev_mps[row], "Link, r") for row in 1:Ny]
    double_hori_A  = [IndexSet(hori_A_inds[row]..., prime(hori_A_inds[row])...) for row in 1:Ny]
    A_prev_unique  = [setdiff(double_hori_A[row], hori_prev_inds[row]) for row in 1:Ny]
    return A_prev_unique
end

function initGuessRandomMPS(A::fPEPS, prev_mps, col::Int, chi::Int)
    @timeit "randomMPS" begin
        Ny, Nx = size(A)
        A_prev_unique  = getGuessInds(A, prev_mps, col)
        hori_cmbs      = Vector{ITensor}(undef, Ny)
        hori_cis       = Vector{Index}(undef, Ny)
        @inbounds for row in 1:Ny
            cmb            = combiner(A_prev_unique[row]..., tags="r$row,CMB,Site")
            hori_cmbs[row] = cmb
            hori_cis[row]  = combinedind(cmb)
        end
        guess = randomMPS(hori_cis, chi)
        @inbounds for row in 1:Ny
            guess[row] *= hori_cmbs[row]
        end
    end
    return guess
end

function initGuessRandomITensor(A::fPEPS, prev_mps, col::Int, chi::Int)
    @timeit "randomITensor" begin
        Ny, Nx = size(A)
        A_prev_unique  = getGuessInds(A, prev_mps, col)
        guess = MPS(Ny)
        up_inds        = [Index(chi, "Link,u,c$col,r$row") for row in 1:Ny-1]
        @inbounds for row in 1:Ny
            row_inds = nothing
            if row == 1
                row_inds = IndexSet(A_prev_unique[row]..., up_inds[row])
            elseif 1 < row < Ny
                row_inds = IndexSet(A_prev_unique[row]..., up_inds[row-1], up_inds[row])
            else
                row_inds = IndexSet(A_prev_unique[row]..., up_inds[row-1])
            end
            #guess[row] = is_cu ? randomCuITensor(row_inds) : randomITensor(row_inds)
            guess[row] = randomITensor(row_inds)
        end
    end
    return guess
end

function constructEnvsAbove(A::fPEPS, guess::MPS, prev_mps::Vector{<:ITensor}, ops::Vector{ITensor}, col::Int, sweep::Int)
    is_cu  = is_gpu(A)
    Ny, Nx = size(A)
    order  = isodd(sweep) ? (1:Ny) : reverse(1:Ny)
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    Envs_above = Vector{ITensor}(undef, Ny+1)
    Envs_above[end]  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @inbounds for (ii, env_row) in enumerate(reverse(order))
        tmp  = A[env_row, col]
        tmp *= ops[env_row]
        tmp *= prev_mps[env_row]
        next_row   = isodd(sweep) ? env_row + 1 : env_row - 1
        Env_above  = next_row > 0 ? copy(Envs_above[next_row]) : dummy
        Env_above *= tmp
        Env_above *= dag(prime(A[env_row, col]))
        Env_above *= guess[env_row]
        Envs_above[env_row] = Env_above
    end
    return Envs_above
end

function fitPEPSMPOold(A::fPEPS, prev_mps::Vector{<:ITensor}, ops::Vector{ITensor}, col::Int, chi::Int)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    @timeit "make guess" begin
        guess::MPS = initGuessRandomMPS(A, prev_mps, col, chi)
        if is_cu
            guess = cuMPS(guess)
        end
    end
    #@show inner(guess, guess)
    dummy = is_cu ? cuITensor(1.0) : ITensor(1.0)
    for sweep in 1:2
        order      = isodd(sweep) ? (1:Ny) : reverse(1:Ny)
        Envs_above = constructEnvsAbove(A, guess, prev_mps, ops, col, sweep)
        Env_below  = is_cu ? cuITensor(1.0) : ITensor(1.0)
        @inbounds for row in order
            # construct environment for the row, have to do this every time
            @timeit "build row Env" begin
                Env  = copy(Env_below)
                tmp  = A[row, col]
                tmp *= ops[row]
                tmp *= prev_mps[row]
                Env *= dag(prime(A[row, col]))
                Env *= tmp
                Env_below  = copy(Env)
                next_row = isodd(sweep) ? row+1 : row -1
                Env *= next_row > 0 ? copy(Envs_above[next_row]) : dummy
                #guess[row] = copy(Env)
                #Env_below  *= guess[row]
            end
            # get l-inds for svd
            @timeit "do SVD" begin
                if isodd(sweep) && row < Ny
                    svd_linds     = setdiff(inds(guess[row]),IndexSet(commonind(guess[row], guess[row+1])))
                    tsvd          = svd(Env, svd_linds, mindim=chi, maxdim=chi)
                    # update guess at row
                    guess[row]    = tsvd.U
                    guess[row+1] *= tsvd.S*tsvd.V
                elseif !isodd(sweep) && row > 1
                    svd_linds     = setdiff(inds(guess[row]),IndexSet(commonind(guess[row], guess[row-1])))
                    tsvd          = svd(Env, svd_linds, mindim=chi, maxdim=chi)
                    # update guess at row
                    guess[row]    = tsvd.U
                    guess[row-1] *= tsvd.S*tsvd.V
                else
                    # update guess at row
                    guess[row]   = copy(Env)
                end
                Env_below    *= guess[row]
            end
        end
    end
    return guess
end

# now with two-site
function fitPEPSMPO(A::fPEPS, prev_mps::Vector{<:ITensor}, ops::Vector{ITensor}, col::Int, chi::Int)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    @timeit "make guess" begin
        guess::MPS = initGuessRandomMPS(A, prev_mps, col, chi)
        if is_cu
            guess = cuMPS(guess)
        end
    end
    for sweep in 1:1
        @timeit "make Envs_above" begin
            Envs_above = constructEnvsAbove(A, guess, prev_mps, ops, col, sweep)
        end
        order     = isodd(sweep) ? (1:Ny-1) : reverse(1:Ny-1)
        Env_below = is_cu ? cuITensor(1.0) : ITensor(1.0)
        @inbounds for row in order
            # construct environment for the row, have to do this every time
            @timeit "build row Env" begin
                Env  = copy(Env_below)
                for row_ in (row, row+1)
                    @timeit "build tmp" begin
                        tmp = A[row_, col]
                        tmp *= ops[row_]
                        tmp *= prev_mps[row_]
                    end
                    @timeit "tmp into Env" begin
                        Env *= tmp
                    end
                    Env *= dag(prime(A[row_, col]))
                    if row_ == row
                        Env_below = copy(Env)
                    end
                end
                Env *= Envs_above[row+2]
            end
            # get l-inds for svd
            @timeit "do SVD" begin
                svd_linds    = setdiff(inds(guess[row]),IndexSet(commonind(guess[row], guess[row+1])))
                tsvd         = svd(Env, svd_linds, mindim=chi, maxdim=chi)
                # update guess at row
                guess[row]   = tsvd.U
                guess[row+1] = tsvd.S*tsvd.V
                Env_below   *= guess[row]
            end
        end
    end
    return guess
end

function buildEdgeEnvironment(A::fPEPS, H, left_H_terms, side::Symbol, col::Int; kwargs...)::Environments
    Ny, Nx          = size(A)
    is_cu::Bool     = is_gpu(A)
    chi::Int        = get(kwargs, :env_maxdim, 1)
    dummy = Vector{ITensor}(undef, Ny)
    for row in 1:Ny
        dummy[row]  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    end
    dummy_mps       = MPS(dummy, 0, Ny+1)
    I_mps           = buildNewI(A, dummy_mps, col, chi)
    field_H_terms   = getDirectional(vcat(H[:, col]...), Field)
    vert_H_terms    = getDirectional(vcat(H[:, col]...), Vertical)
    vHs             = [buildNewVerticals(A, vert_H_terms[vert_op], dummy_mps, col, chi) for vert_op in 1:length(vert_H_terms)]
    fHs             = [buildNewFields(A, field_H_terms[field_op], dummy_mps, col, chi) for field_op in 1:length(field_H_terms)]
    Hs              = vcat(vHs, fHs)
    maxdim::Int     = get(kwargs, :maxdim, 1)
    cutoff::Float64 = get(kwargs, :cutoff, 0.0)
    hori_inds = Vector{IndexSet}(undef, Ny)
    hori_cmbs = Vector{ITensor}(undef, Ny)
    up_inds   = [Index(chi, "Link,c$col,r$row,u") for row in 1:Ny-1]
    @inbounds for row in 1:Ny
        hori_inds[row] = intersect(map(x->inds(x[row], "Link,r"), Hs))[1]
        hori_cmbs[row] = combiner(hori_inds[row], tags="r$row,CMB,Site")
        for ii in 1:length(Hs)
            Hs[ii][row] *= hori_cmbs[row]
        end
    end
    H_overall       = sum(Hs; cutoff=cutoff, maxdim=chi)
    @inbounds for row in 1:Ny
        H_overall[row] *= hori_cmbs[row]
    end
    side_H          = side == :left ? H[:, col] : H[:, col - 1]
    side_H_terms    = getDirectional(vcat(side_H...), Horizontal)
    in_progress     = Matrix{ITensor}(undef, Ny, length(side_H_terms))
    @inbounds for side_term in 1:length(side_H_terms)
        in_progress[1:Ny, side_term] = generateEdgeDanglingBonds(A, side_H_terms[side_term], side, col, chi)
    end
    diag_H_terms     = getDirectional(vcat(side_H...), Diag)
    diag_in_progress = Matrix{ITensor}(undef, Ny, length(diag_H_terms))
    @inbounds for diag_term in 1:length(diag_H_terms)
        diag_in_progress[1:Ny, diag_term] = generateEdgeDanglingBonds(A, diag_H_terms[diag_term], side, col, chi)
    end
    return Environments(I_mps, H_overall, in_progress, diag_in_progress)
end

function buildNextEnvironment(A::fPEPS, prev_Env::Environments, H,
                              side::Symbol,
                              col::Int;
                              kwargs...)
    Ny, Nx        = size(A)
    chi::Int      = get(kwargs, :env_maxdim, 1)
    @timeit "build new_I and new_H" begin
        new_I = buildNewI(A, prev_Env.I, col, chi)
        new_H = buildNewI(A, prev_Env.H, col, chi)
    end
    field_H_terms = getDirectional(vcat(H[:, col]...), Field)
    vert_H_terms  = getDirectional(vcat(H[:, col]...), Vertical)
    hori_H_terms  = getDirectional(vcat(H[:, col]...), Horizontal)
    diag_H_terms  = getDirectional(vcat(H[:, col]...), Diag)
    side_H        = side == :left ? H[:, col] : H[:, col - 1]
    side_H_terms  = getDirectional(vcat(side_H...), Horizontal)
    diag_side_H_terms = getDirectional(vcat(side_H...), Diag)
    H_term_count  = 1 + length(field_H_terms) + length(vert_H_terms)
    H_term_count += (side == :left ? length(side_H_terms) : length(hori_H_terms))
    H_term_count += (side == :left ? length(diag_side_H_terms) : length(diag_H_terms))
    new_H_mps     = Vector{MPS}(undef, H_term_count)
    new_H_mps[1]  = deepcopy(new_H)
    @timeit "build new verts" begin
        new_H_mps[2:1+length(vert_H_terms)] = [buildNewVerticals(A, vert_H_terms[vert_op], prev_Env.I, col, chi) for vert_op in 1:length(vert_H_terms)]
    end
    @timeit "build new fields" begin
        new_H_mps[2+length(vert_H_terms):1+length(vert_H_terms)+length(field_H_terms)]  = [buildNewFields(A, field_H_terms[field_op], prev_Env.I, col, chi) for field_op in 1:length(field_H_terms)]
    end
    connect_H    = side == :left ? side_H_terms : hori_H_terms
    @timeit "connect dangling bonds" begin
        @inbounds for (cc, cH) in enumerate(connect_H)
            new_H = connectDanglingBonds(A, cH, prev_Env.InProgress[:, cc], side, col; kwargs...)
            new_H_mps[length(vert_H_terms) + length(field_H_terms) + 1 + cc] = MPS(new_H, 0, Ny+1)
        end
    end
    diag_connect_H = side == :left ? diag_side_H_terms : diag_H_terms
    @timeit "connect dangling bonds" begin
        @inbounds for (cc, cH) in enumerate(diag_connect_H)
            new_H = connectDanglingBonds(A, cH, prev_Env.DiagInProgress[:, cc], side, col; kwargs...)
            new_H_mps[length(vert_H_terms) + length(field_H_terms) + length(connect_H) + 1 + cc] = MPS(new_H, 0, Ny+1)
        end
    end
    maxdim::Int     = get(kwargs, :maxdim, 1)
    cutoff::Float64 = get(kwargs, :cutoff, 0.0)
    @timeit "sum H mps" begin
        hori_inds = Vector{IndexSet}(undef, Ny)
        hori_cmbs = Vector{ITensor}(undef, Ny)
        up_inds   = [Index(chi, "Link,c$col,r$row,u") for row in 1:Ny-1]
        for row in 1:Ny
            hori_inds[row] = intersect(map(x->inds(x[row], "Link,r"), new_H_mps))[1]
            hori_cmbs[row] = combiner(hori_inds[row], tags="r$row,CMB,Site")
            for ii in 1:length(new_H_mps)
                new_H_mps[ii][row] *= hori_cmbs[row]
            end
        end
        H_overall    = sum(new_H_mps; cutoff=cutoff, maxdim=chi)
        for row in 1:Ny
            H_overall[row] *= hori_cmbs[row]
        end
    end
    gen_H_terms  = side == :left ? hori_H_terms : side_H_terms
    in_progress  = Matrix{ITensor}(undef, Ny, length(gen_H_terms))
    @timeit "gen dangling bonds" begin
        @inbounds for side_term in 1:length(gen_H_terms)
            in_progress[1:Ny, side_term] = generateNextDanglingBonds(A, gen_H_terms[side_term], prev_Env.I, side, col; kwargs...)
        end
    end
    diag_gen_H_terms  = side == :left ? diag_H_terms : diag_side_H_terms
    diag_in_progress  = Matrix{ITensor}(undef, Ny, length(diag_gen_H_terms))
    @timeit "diag gen dangling bonds" begin
        @inbounds for side_term in 1:length(diag_gen_H_terms)
            diag_in_progress[1:Ny, side_term] = generateNextDanglingBonds(A, diag_gen_H_terms[side_term], prev_Env.I, side, col; kwargs...)
        end
    end
    return Environments(new_I, H_overall, in_progress, diag_in_progress)
end

function buildNewVerticals(A::fPEPS, H, prev_I::MPS, col::Int, chi::Int)::MPS
    Ny, Nx              = size(A)
    is_cu               = is_gpu(A)
    col_site_inds       = [firstind(A[row, col], "Site") for row in 1:Ny]
    ops                 = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    vertical_row_a      = H.sites[1][1]
    vertical_row_b      = H.sites[2][1]
    ops[vertical_row_a] = replaceind!(copy(H.ops[1]), H.site_ind, col_site_inds[vertical_row_a])
    ops[vertical_row_a] = replaceind!(ops[vertical_row_a], H.site_ind', col_site_inds[vertical_row_a]')
    ops[vertical_row_b] = replaceind!(copy(H.ops[2]), H.site_ind, col_site_inds[vertical_row_b])
    ops[vertical_row_b] = replaceind!(ops[vertical_row_b], H.site_ind', col_site_inds[vertical_row_b]')
    #return fitPEPSMPOold(A, data(prev_I), ops, col, chi)
    return fitPEPSMPO(A, data(prev_I), ops, col, chi)
end

function buildNewFields(A::fPEPS, H, prev_I::MPS, col::Int, chi::Int)::MPS
    Ny, Nx         = size(A)
    is_cu          = is_gpu(A)
    col_site_inds  = [firstind(A[row, col], "Site") for row in 1:Ny]
    ops            = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    field_row      = H.sites[1][1]
    ops[field_row] = replaceind!(copy(H.ops[1]), H.site_ind, col_site_inds[field_row])
    ops[field_row] = replaceind!(ops[field_row], H.site_ind', col_site_inds[field_row]')
    #return fitPEPSMPOold(A, data(prev_I), ops, col, chi)
    return fitPEPSMPO(A, data(prev_I), ops, col, chi)
end

function buildNewI(A::fPEPS, prev_I::MPS, col::Int, chi::Int)::MPS
    Ny, Nx         = size(A)
    is_cu          = is_gpu(A)
    col_site_inds  = [firstind(A[row, col], "Site") for row in 1:Ny]
    ops            = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    #return fitPEPSMPOold(A, data(prev_I), ops, col, chi)
    return fitPEPSMPO(A, data(prev_I), ops, col, chi)
end

function generateNextDanglingBonds(A::fPEPS,
                                   H,
                                   Ident::MPS,
                                   side::Symbol,
                                   col::Int;
                                   kwargs...)::Vector{ITensor}
    Ny, Nx          = size(A)
    is_cu           = is_gpu(A)
    chi::Int        = get(kwargs, :env_maxdim, 1)
    op_row          = side == :left ? H.sites[1][1] : H.sites[2][1]
    H_op            = side == :left ? H.ops[1]      : H.ops[2]
    col_site_inds   = [firstind(A[row, col], "Site") for row in 1:Ny]
    ops             = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    ops[op_row]     = replaceind!(copy(H_op), H.site_ind, col_site_inds[op_row])
    ops[op_row]     = replaceind!(ops[op_row], H.site_ind', col_site_inds[op_row]')
    #return data(fitPEPSMPOold(A, data(Ident), ops, col, chi))
    return data(fitPEPSMPO(A, data(Ident), ops, col, chi))
end

function generateEdgeDanglingBonds(A::fPEPS, H, side::Symbol, col::Int, chi::Int)::Vector{ITensor}
    Ny, Nx         = size(A)
    is_cu          = is_gpu(A)
    dummy          = [is_cu ? cuITensor(1.0) : ITensor(1.0) for row in 1:Ny]
    op_row         = side == :left ? H.sites[1][1] : H.sites[2][1]
    H_op           = side == :left ? H.ops[1]      : H.ops[2]
    col_site_inds  = [firstind(A[row, col], "Site") for row in 1:Ny]
    ops            = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    ops[op_row]    = replaceind!(copy(H_op), H.site_ind, col_site_inds[op_row])
    ops[op_row]    = replaceind!(ops[op_row], H.site_ind', col_site_inds[op_row]')
    #return data(fitPEPSMPOold(A, dummy, ops, col, chi))
    return data(fitPEPSMPO(A, dummy, ops, col, chi))
end

function connectDanglingBonds(A::fPEPS,
                              oldH,
                              in_progress::Vector{ITensor},
                              side::Symbol,
                              col::Int;
                              kwargs...)::Vector{ITensor}
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    chi::Int = get(kwargs, :env_maxdim, 1)
    op_row_a = oldH.sites[1][1]
    op_row_b = oldH.sites[2][1]
    op       = side == :left ? oldH.ops[2] : oldH.ops[1]
    application_row      = side == :left ? op_row_b : op_row_a
    col_site_inds        = [firstind(A[row, col], "Site") for row in 1:Ny]
    ops                  = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    ops[application_row] = replaceind!(copy(op), oldH.site_ind, col_site_inds[application_row])
    ops[application_row] = replaceind!(ops[application_row], oldH.site_ind', col_site_inds[application_row]')
    in_prog_mps          = MPS(in_progress, 0, Ny + 1)
    #return data(fitPEPSMPOold(A, in_progress, ops, col, chi))
    return data(fitPEPSMPO(A, in_progress, ops, col, chi))
end

function buildLs(A::fPEPS, H; kwargs...)
    Ny, Nx         = size(A)
    Ls             = Vector{Environments}(undef, Nx)
    start_col::Int = get(kwargs, :start_col, 1)
    if start_col == 1
        left_H_terms = getDirectional(vcat(H[:, 1]...), Horizontal)
        Ls[1]        = buildEdgeEnvironment(A, H, left_H_terms, :left, 1; kwargs...)
    end
    loop_col = start_col == 1 ? 2 : start_col
    @inbounds for col in loop_col:(Nx-1)
        @debug "Building left col $col"
        Ls[col] = buildNextEnvironment(A, Ls[col-1], H, :left, col; kwargs...)
    end
    return Ls
end

function buildRs(A::fPEPS, H; kwargs...)
    Ny, Nx         = size(A)
    start_col::Int = get(kwargs, :start_col, Nx)
    Rs             = Vector{Environments}(undef, Nx)
    if start_col == Nx
        right_H_terms = getDirectional(vcat(H[:, Nx-1]...), Horizontal)
        Rs[Nx]        = buildEdgeEnvironment(A, H, right_H_terms, :right, Nx; kwargs...)
    end
    loop_col = start_col == Nx ? Nx - 1 : start_col
    @inbounds for col in reverse(2:loop_col)
        @debug "Building right col $col"
        Rs[col] = buildNextEnvironment(A, Rs[col+1], H, :right, col; kwargs...)
    end
    return Rs
end


text:
using econ627ps7
using Test
g(x)=(x-1.0)^2
gprime(x)=2*(x-1)

@testset "Basic Newtons Method Tests" begin
    @test newton(g,gprime,10.0) ≈ 1.0  atol=0.00001
    @test newton(g,10)≈ 1.0 atol=0.00001
end;

@testset "Convergence" begin
    @test (newton(x->x^2+1,x₀=10)==nothing)
    @test (newton(x->x^2+1,x->2*x,x₀=10)==nothing)
end;

@testset "Max Iter" begin
    @test (newton(g,x₀=10,maxiter=5)==nothing)
    @test (newton(g,gprime,x₀=10,maxiter=5)==nothing)
end;

@testset "Tolerance" begin
    x1=(newton(g,x₀=10,tol=1E-6))
    x2=(newton(g,x₀=10,tol=1E-5))
    @test(abs(x1-1)<abs(x2-1))
end;

@testset "BigFloat" begin
    @test (newton(g,gprime,x₀=BigFloat(10))≈1.0)
    @test (newton(g,x₀=BigFloat(10))≈1.0)
end;


text:
@testset "fsm_stream_reassembler_cap.jl" begin
    @testset begin
        rs = StreamReassembler(2)
        push_substring!(rs, "ab", 0)
        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"

        push_substring!(rs, "cd", 2)
        @test rs.cur_idx == 4
        @test read!(rs.output) == "cd"

        push_substring!(rs, "ef", 4)
        @test rs.cur_idx == 6
        @test read!(rs.output) == "ef"
    end

    @testset begin
        rs = StreamReassembler(2)
        push_substring!(rs, "ab", 0)
        @test rs.cur_idx == 2
        push_substring!(rs, "cd", 2)

        @test rs.cur_idx == 2
        @test read!(rs.output) == "ab"

        push_substring!(rs, "cd", 2)
        @test rs.cur_idx == 4
        @test read!(rs.output) == "cd"
    end
    @testset begin
        rs = StreamReassembler(1)
        push_substring!(rs, "ab", 0)
        @test rs.cur_idx == 1

        push_substring!(rs, "ab", 0)
        @test rs.cur_idx == 1

        @test read!(rs.output) == "a"

        push_substring!(rs, "abc", 0)
        @test rs.cur_idx == 2
        @test read!(rs.output) == "b"
        @test rs.cur_idx == 2
    end

    @testset begin
        rs = StreamReassembler(3)
        for i in 0:3:99997
            seg = Char.(rand(1:125, 6)) |> string
            push_substring!(rs, seg, i)
            @test rs.cur_idx == i + 3
            @test read!(rs.output) == SubString(seg, 1:3)
        end
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
module Threshold

using CuArrays

include("poisson.jl")

end


text:
#=
Multi-armed Bandit problem

Simulates the k-armed bandit problem
=#

using Distributions
using PyPlot
using Statistics

function karm_bandit(k::Int64=10, steps::Int64=1000, runs::Int64=1)
# k : Number of arms
    # mean_range : range of values that the mean for each arm is sampled from
    # steps : number of steps the simulation runs for.

    # Array of distributions
    means = [m for m in rand(Normal(),k)]
    distributions = [Normal(m) for m in means]

    greedy_rewards = Array{Float64}(undef,steps,runs)
    eps1_rewards = Array{Float64}(undef,steps,runs)
    eps2_rewards = Array{Float64}(undef,steps,runs)
    # cumulative_rewards = Array{Float64}(undef,steps)
    # cumulative_average = Array{Float64}(undef,steps)

    eps1 = .1
    eps2 = .01

    # Main loop
    for jj=1:runs
        # Q array of reward estimates
        Q_greedy = [0.0 for i in 1:k]
        Q_eps1 = [0.0 for i in 1:k]
        Q_eps2 = [0.0 for i in 1:k]

        # Counts how many times Q has been called (n):
        Qn_greedy = [0 for i in 1:k]
        Qn_eps1 = [0 for i in 1:k]
        Qn_eps2 = [0 for i in 1:k]

        # Simulate rewards
        # ith row is ith time step, jth column is jth run

        for ii=1:steps

            # Choose action
            action_greedy = totally_greedy(Q_greedy,Qn_greedy)
            action_eps1 = explore(Q_eps1,Qn_eps1,eps1)
            action_eps2 = explore(Q_eps2,Qn_eps2,eps2)

            # Store rewards
            greedy_rewards[ii,jj] = rand(distributions[action_greedy])
            eps1_rewards[ii,jj] = rand(distributions[action_eps1])
            eps2_rewards[ii,jj] = rand(distributions[action_eps2])

            # println(rewards[ii,jj])

            # Update estimate of q*
            if Qn_greedy[action_greedy] > 0
                # Normal update
                Q_greedy[action_greedy] += 1/Qn_greedy[action_greedy]*(greedy_rewards[ii,jj] - Q_greedy[action_greedy])
                Qn_greedy[action_greedy] += 1
            else
                # Just increment Qn_eps1
                Q_greedy[action_greedy] = greedy_rewards[ii,jj]
                Qn_greedy[action_greedy] += 1
            end

            if Qn_eps1[action_eps1] > 0
                # Normal update
                Q_eps1[action_eps1] += 1/Qn_eps1[action_eps1]*(eps1_rewards[ii,jj] - Q_eps1[action_eps1])
                Qn_eps1[action_eps1] += 1
            else
                # Just increment Qn
                Q_eps1[action_eps1] = eps1_rewards[ii,jj]
                Qn_eps1[action_eps1] += 1
            end

            if Qn_eps2[action_eps2] > 0
                # Normal update
                Q_eps2[action_eps2] += 1/Qn_eps2[action_eps2]*(eps2_rewards[ii,jj] - Q_eps2[action_eps2])
                Qn_eps2[action_eps2] += 1
            else
                # Just increment Qn
                Q_eps2[action_eps2] = eps2_rewards[ii,jj]
                Qn_eps2[action_eps2] += 1
            end

            # if ii == 1
            #     cumulative_rewards[ii] = rewards[ii]
            #     cumulative_average[ii] = rewards[ii]
            # else
            #     cumulative_rewards[ii] = cumulative_rewards[ii-1] + rewards[ii]
            #     cumulative_average[ii] = cumulative_rewards[ii] / ii
            # end
        end
    end

    # Create average reward over all runs
    avg_rewards_greedy = Array{Float64}(undef,steps)
    avg_rewards_eps1 = Array{Float64}(undef,steps)
    avg_rewards_eps2 = Array{Float64}(undef,steps)

    for kk=1:steps
        avg_rewards_greedy[kk] = mean(greedy_rewards[kk,:])
        avg_rewards_eps1[kk] = mean(eps1_rewards[kk,:])
        avg_rewards_eps2[kk] = mean(eps2_rewards[kk,:])
    end

    plot(1:(steps+1), vcat(0.0,avg_rewards_greedy), linewidth=1, label="Greedy");
    plot(1:(steps+1), vcat(0.0,avg_rewards_eps1), linewidth=1, label=string("Eps1 = ", eps1));
    plot(1:(steps+1), vcat(0.0,avg_rewards_eps2), linewidth=1, label=string("Eps2 = ", eps2));

    # Max reward
    plot(1:(steps+1), [maximum(means) for i in 1:(steps+1)], linewidth=1, linestyle="--", label="Maximum")
    # display(sum_plot)

    legend()
    means

end

function totally_greedy(Q,Qn)::Int64
    # Make sure all options are chosen at least once
    # println(length(findall(Qn .== 0)))
    # if length(findall(Qn .== 0)) > 0
    #     unexplored_vec = findall(Qn .== 0)
    #     # println("length(findall(Qn .== 0)) < 0")
    #     return unexplored_vec[rand(1:length(unexplored_vec))]
    # else
        if length(findall(Q .== maximum(Q))) == 1
            # println("length(findall(Q .== maximum(Q))) == 1")
            return argmax(Q)
        else
            # Choose randomly among largest values
            # println("else")
            max_vec = findall(Q .== maximum(Q))
            return max_vec[rand(1:length(max_vec))]
        end
    # end
end

function explore(Q,Qn,eps)::Int64
    explore_bool = rand(Binomial(1,eps)) # Takes exploration action with probability 1 - eps. Binomial
    if explore_bool == 1
        # Exploration action
        # println("explore_bool == 1")
        return rand(1:length(Q))
    else
        # Greedy action
        # println("greedy action")
        # println("START")
        # println(Q)
        max_vec = findall(Q .== maximum(Q))
        # println(max_vec)
        # println(max_vec[rand(1:length(max_vec))])
        # println("END")
        return max_vec[rand(1:length(max_vec))]
    end
end


text:
# SPDX-License-Identifier: X11
# 2020-09-13
# Dist Max (500pt)

function _main()
  N = parse(Int, readline())
  X, Y = parse.(Int, split(readline()))
  A = [(X + Y, X - Y)]
  for i ∈ 2:N
    X, Y = parse.(Int, split(readline()))
    push!(A, (X + Y, X - Y))
  end
  return max(maximum(x -> x[1], A) - minimum(x -> x[1], A),
             maximum(x -> x[2], A) - minimum(x -> x[2], A))
end

println(_main())

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module TriMeshRender

using TriMeshes
using Luxor

export isometric, render,  box!

quad! = TriMeshes.quad!

const r6 = sqrt(6)
const rr6 = 1/r6
const r2 = sqrt(2)
const r3 = sqrt(3)

const rM = [ [r3, 1, r2] [0,2,r2] [-r3,1,r2] ]
const Orth = [ [1,0,0] [0,1,0] [0,0,0] ]

function isometric(v::Vertex)
	iv = Orth * rr6 * rM * v
	Point(iv[1], iv[2])
end

isometric(n::Net) = map(isometric, n.vertices)

function render(n::Net, w, h, fname)
	@png begin
	points = isometric(n)
	for f in n.faces
		line(points[f.AB.from], points[f.AB.to], :stroke)
		line(points[f.BC.from], points[f.BC.to], :stroke)
		line(points[f.CA.from], points[f.CA.to], :stroke)
	end
	end w h fname
end

###
end


text:
export wignersurmise,unfoldedeigvals ,unfoldedSpacings

"""
```julia
wignersurmise(s ; beta)
```
- `beta` : default `1`,

    `1` for the spacing density corresponds to GOE(2)

    ``p_1(s)=\\frac{\\pi s}{2} e^{-\\pi s^{2} / 4}``

    `2` for the spacing density corresponds to GUE(2)

    ``p_2(s)=\\frac{32 s^{2}}{\\pi^{2}} e^{-4 s^{2} / \\pi}``

-  For more info see [Wigner surmise on Wikipedia](https://en.wikipedia.org/wiki/Wigner_surmise)

# Examples
Evaluate Wigner surmise for GOE at `2`
```julia
wignersurmise(2)

0.1357605281502967
```

Evaluate Wigner surmise for GUE at `2`
```julia
wignersurmise(2,beta=2)

0.07962814352725955
```
"""
function wignersurmise(s::Real; beta = 1 :: Int)
    if beta == 1
        return (pi/2) * s * exp(-pi * s^2 / 4)
    elseif beta == 2
        return 32*(s/pi)^2 * exp(- 4 * s^2 / pi)
    end
end

"""
```julia
unfoldedeigvals(E; deg, goe)
```
- `E` : The array of eigenvalues
- `deg` : default `10`
- `goe` : default `false`
"""
function unfoldedeigvals(E::Array; deg=10::Int, goe = false::Bool)
    if goe
        return length(E)*cdf(Semicircle(2),E)
    end

    return polyfit(E,(x->sum(E.<x)).(E),deg).(E)

end
function unfoldedSpacings(E::Array;deg=10::Int, goe = false::Bool)
    return diff(unfoldedeigvals(E,deg=deg,goe=goe))
end


text:
module Preconditioners

using AlgebraicMultigrid

using JuLIP: AbstractAtoms, Dofs, maxdist,
            cutoff, positions, forces, mat, vecs,
            set_positions!, chemical_symbols, rnn, JVec, JMat ,
            AbstractCalculator, calculator, cell,
            fixedcell, projectxfree, _pos_to_dof

using JuLIP.Potentials: @pot, @analytic, evaluate, evaluate_d, PairPotential, HS,
                        SitePotential, sites, C0Shift, _precon_or_hessian_pos,
                        MSitePotential

using SparseArrays: SparseMatrixCSC

using LinearAlgebra: cholesky, I, Symmetric, norm

import SuiteSparse

import JuLIP:             update!
import JuLIP.Potentials:  precon!, cutoff
import Base:              *, \, size
import LinearAlgebra:     ldiv!, mul!, dot

export Exp, FF

"""
`IPPrecon`: the standard preconditioner type for JuLIP

`IPPrecon` stores a field `p <: AbstractCalculator` which is
should implement `precon(p, r, R)` which is a block matrix analogue of
`hess(p, r, R)` and is used to assemble the preconditioner matrix via
```julia
precon_matrix(P, at::AbstractAtoms)
```
If `p` should be updated in response to an update of `at` then
the type can overload `JuLIP.Preconditioners.update_inner!`.

## Constructor:
```
IPPrecon(p::AbstractCalculator, at::AbstractAtoms; kwargs...)
```

## Keyword arguments:

* `updatedist`: determines after how much movement of the atoms, the
preconditioner is updated
* `tol`: solver tolerance if AMG is used
* `updatefreq`: after how many iterations do we force an update {10}
* `stab`: stabilisation constant {0.01}, add `stab * I` to `P`
* `solve`: which solver to used, `:amg` or `:chol`, default is `:amg`
"""
mutable struct IPPrecon{TV, TS, T <: AbstractFloat, TI <: Integer}
   p::TV
   solver::TS
   A::SparseMatrixCSC{T, TI}
   oldX::Vector{JVec{T}}
   updatedist::T
   tol::T
   updatefreq::TI
   skippedupdates::TI
   stab::T
   innerstab::T
end

function IPPrecon(p, at::AbstractAtoms;
         updatedist=0.2 * rnn(at), tol=1e-7, updatefreq=10, stab=0.01,
         solver = :chol, innerstab=0.0)
   # make sure we don't use this in a context it is not intended for!
   @assert fixedcell(at)
   A = AlgebraicMultigrid.poisson(12)
   if solver == :amg
      solver = ruge_stuben(A)
   elseif solver in [:chol, :direct]
      solver = cholesky(A)
   else
      error("`IPPrecon` : unknown solver $(solver)")
   end
   P = IPPrecon(p, solver, A, positions(at), updatedist, tol, updatefreq, 0, stab, innerstab)
   # the force_update! makes the first assembly pass.
   return force_update!(P, at)
end

# some standard Base functionality lifted to IPPrecon
ldiv!(out::Dofs, P::IPPrecon{TV, TS}, f::Dofs) where TV where TS <: AlgebraicMultigrid.MultiLevel =
   copyto!(out, solve(P.solver, f, AlgebraicMultigrid.V(); maxiter = 200, tol = P.tol))
ldiv!(out::Dofs, P::IPPrecon{TV, TS}, f::Dofs) where TV where TS <: SuiteSparse.CHOLMOD.Factor =
   copyto!(out, P.solver \ f)
mul!(out::Dofs, P::IPPrecon, x::Dofs) = mul!(out, P.A, x)
dot(x, P::IPPrecon, y) = dot(x, P * y)
*(P::IPPrecon, x::AbstractVector) = P.A * x
\(P::IPPrecon, x::AbstractVector) = ldiv!(zeros(length(x)), P, x)
Base.size(P::IPPrecon) = size(P.A)

need_update(P::IPPrecon, at::AbstractAtoms) =
   (P.skippedupdates > P.updatefreq) ||
   (maxdist(positions(at), P.oldX) >= P.updatedist)

update!(P::IPPrecon, at::AbstractAtoms) =
   need_update(P, at) ? force_update!(P, at) : (P.skippedupdates += 1; P)

update_solver!(P::IPPrecon{TV, TS}) where TV where TS <: SuiteSparse.CHOLMOD.Factor =
   cholesky(Symmetric(P.A))

update_solver!(P::IPPrecon{TV, TS}) where TV where TS <: AlgebraicMultigrid.MultiLevel =
   ruge_stuben(P.A)

function force_update!(P::IPPrecon, at::AbstractAtoms)
   # perform updates of the potential p (if needed; usually not)
   P.p = update_inner!(P.p, at)
   # construct the preconditioner matrix ...
   Pmat = precon_matrix(P.p, at, P.innerstab)
   Pmat = Pmat + P.stab * I
   Pmat = projectxfree(at, Pmat)
   # and the AMG solver
   P.A = Pmat
   P.solver = update_solver!(P)
   # remember the atom positions
   copyto!(P.oldX, positions(at))
   # and remember that we just did a full update
   P.skippedupdates = 0
   return P
end

# ============== Implementation of some concrete types ======================

# default implementation of update_inner!;
# most of the time no inner update is required

"update `IPPrecon.p`, must return updated `p`"
update_inner!(p::Any, at::AbstractAtoms) = p

"return the three linear indices associated with an atom index"
atind2lininds(i::Integer) = (i-1) * 3 + [1;2;3]

"""
build the preconditioner matrix associated with the potential V
"""
precon_matrix(V, at::AbstractAtoms, innerstab = 0.1;
              preconmap = (hEs, tmp, V, R) -> precon!(hEs, tmp, V, R, innerstab)) =
   _pos_to_dof(_precon_or_hessian_pos(V, at, preconmap), at)

precon_matrix(V::MSitePotential, at::AbstractAtoms, innerstab = 0.1;
              preconmap = (hEs, tmp, V, R, Z, z0) -> precon!(hEs, tmp, V, R, Z, z0, innerstab)) =
   _pos_to_dof(_precon_or_hessian_pos(V, at, preconmap), at)

"""
A variant of the `Exp` preconditioner; see

### Constructor: `Exp(at::AbstractAtoms)`

Keyword arguments:

* `A`: stiffness of potential {default 3.0}
* `r0`: nn distance estimate {default is an automicatic estimate}
* `cutoff_mult`: cut-off is `r0 * cutoff_mult` {default 2.2}
* `kwargs...`: `IPPrecon` parameters

### Reference

      D. Packwood, J. Kermode, L. Mones, N. Bernstein, J. Woolley, N. I. M. Gould,
      C. Ortner, and G. Csanyi. A universal preconditioner for simulating condensed
      phase materials. J. Chem. Phys., 144, 2016.
"""
mutable struct Exp{T, TV}
   Vexp::TV
   energyscale::T
end

cutoff(P::Exp) = cutoff(P.Vexp)

function precon!(hEs, tmp, P::Exp{T}, R::AbstractVector{<: JVec}, innerstab=0) where {T}
   n = length(R)
   for i = 1:n
      hEs[i,i] = (P.energyscale * P.Vexp(norm(R[i]))+innerstab) * one(JMat{T})
   end
   return hEs
end

function Exp(at::AbstractAtoms{T};
             A=T(3.0), r0=rnn(at), cutoff_mult=T(2.2), energyscale = T(1.0),
             kwargs...) where {T}
   e0 = energyscale == :auto ? T(1.0) : energyscale
   rcut = r0 * cutoff_mult
   Vexp = let A=A, r0=r0
      @analytic r -> exp( - A * (r/r0 - 1))
   end * C0Shift(rcut)
   P = IPPrecon(Exp(Vexp, e0), at; kwargs...)
   if energyscale == :auto
      P.p.energyscale = estimate_energyscale(at, P)
      force_update!(P, at)
   end
   return P
end

# want μ * <P v, v> ~ <∇E(x+hv) - ∇E(x), v> / h
function estimate_energyscale(at::AbstractAtoms{T}, P) where {T}
   # get the P-matrix at current configuration
   A = precon_matrix(P.p, at)
   # determine direction in which the cell is maximal
   F = Matrix(cell(at)')
   X0 = positions(at)
   _, i = maximum( (norm(F[:, j]), j) for j = 1:3 ) # hack to make findmax work again
   Fi = JVec{T}(F[:, i])
   # an associated perturbation
   V = [sin(2*pi * dot(Fi, x)/dot(Fi,Fi)) * Fi for x in X0]
   # compute gradient at current positions
   f0 = forces(at) |> mat
   # compute gradient at perturbed positions
   h = 1e-3
   set_positions!(at, X0 + h * V)
   f1 = forces(at) |> mat
   # return the estimated value for the energyscale
   V = mat(V)
   μ = - dot((f1 - f0)/h, V) / dot(A * V[:], V[:])
   if μ < 0.1 || μ > 100.0
      warn("""
      e0 = $(μ) in `estimate_energyscale`; this is likely due to a poor
      starting guess in an optimisation, probably best to set the
      energyscale manually, or use CG instead of LBFGS.
      """)
   end
   return μ
end

"""
`FF`: defines a preconditioner based on a force-field; implementation
is close to the paper where this idea is described:

   Preconditioners for the geometry optimisation and saddle point search of
   molecular systems; Letif Mones, Christoph Ortner & Gábor Csányi;
   Scientific Reports 8, Article number: 13991 (2018)

Each potential has to define a `precon` method from which the preconditioner
is build.
"""
FF(at::AbstractAtoms, V::AbstractCalculator; kwargs...) =
      IPPrecon(V, at; kwargs...)

FF(at::AbstractAtoms; kwargs...) =
      FF(at, calculator(at); kwargs...)

end # end module Preconditioners


text:
#=
Iteration in Julia 0.7
In Julia 0.7, the iteration interface is now just one function: `iterate`.
=#

iterable = [1,2,3,4,5]
@show typeof(iterable)

for a in iterable
    @show a
end

# this is equivalent as the expressin using `while` loop:

iter_result = iterate(iterable)
while iter_result !== nothing
    global iter_result
    (element, state) = iter_result
    # body
    @show element state
    iter_result = iterate(iterable, state)
end

struct EveryNthStep
    start::Int
    length::Int
    step::Int
    endpoint::Bool
    EveryNthStep(;start,length,step) = new(start,length,step)
end

function Base.iterate(iter::EveryNthStep, state=(iter.start, 0))
    element, count = state
    @show state
    if element >= iter.length
        return nothing
    end

    return (element, (element + iter.step, count + 1))
end

#function Base.iterate(it::EveryNth, (element, count)=(it.start, 0))
#    return count >= it.length ? nothing : (element, (element + it.n, count + 1))
#end
Base.eltype(iter::EveryNthStep) = Int

for element in EveryNthStep(start=0,length=10,step=2)
    println(element)
end

for element in EveryNthStep(-30, 30, 2)
    println(element)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using GPX
using GPX: read_gpx_file

fname = joinpath("..", "sample", "data", "sample.gpx")
println(fname)
itr = read_gpx_file(fname)
# positions = collect(itr)
# println(positions)


text:
# arules_testing.jl

using RCall
using DataFrames
using StatsBase                        # for sample() function

include("arules.jl")
include("utils.jl")

init_arules()                       # this is called automatically by `using AssociationRules`

# simulate transactions
groceries = ["milk", "bread", "eggs", "apples", "oranges", "beer"]
transacts = [sample(groceries, 4, replace = false) for x in 1:1000]
transacts_df = DataFrame(transactions(transacts))

# minimum support of 0.2, minimum confidence of 0.8
rules = apriori(transacts_df, 0.2, 0.8)

adult = dataset("adult")            # example dataset included in package
sanitize_input!(adult)              # remove excluded character patterns (e.g., " => ")

# minimum support of 0.1, minimum confidence of 0.8
rules = apriori(adult, 0.1, 0.8)


text:
module plotfxn

    using GeometryTypes
    using Plots
    using algfxn
    using algT

    ###  Plots.jl recipes

    @recipe function f(r::HyperRectangle)
        points = decompose(Point{2,Float64}, r)
        rectpoints = points[[1,2,4,3],:]
        xs = [pt[1] for pt in rectpoints];
        ys = [pt[2] for pt in rectpoints];
        seriestype := :shape
        color := :orange
        linecolor := :black
        #m = (:black, stroke(0))
        s = Shape(xs[:], ys[:])
    end

    @recipe function f(pt::Point)
        xs = [pt[1]]
        ys = [pt[2]]
        seriestype --> :scatter
        color := :orange
        markersize := 3
        xs, ys
    end

    @recipe function f(l::LineSegment)
        xs = [ l[1][1], l[2][1] ]
        ys = [ l[1][2], l[2][2] ]
        # seriestype = :line
        color := :darkblue
        lw := 3
        xs, ys
    end

    @recipe function f(rectList::Vector{<:HyperRectangle})
        for r in rectList
            @series begin
                r
            end
        end
    end

    @recipe function f(ptList::Vector{<:Point})
        for p in ptList
            @series begin
                p
            end
        end
    end

    @recipe function f(lineList::Vector{<:LineSegment})
        for l in lineList
            @series begin
                l
            end
        end

    end

    function plotRoom(room)
        aPlot = plot() #Todo! this assumes plotroom is first thing called()
        roomWidth, roomHeight, walls, obstacles = room.width, room.height, room.walls, room.obstacles
        plot!(aPlot, walls, color =:black)
        plot!(aPlot, obstacles, fillalpha=0.5)
        return aPlot
    end

    function plotRRT(roomPlot, graph, solPath, title)
        startstate, goalstate, nodeslist, edgeslist = graph.startstate, graph.goalstate, graph.nodeslist, graph.edgeslist

        x = [n.state[1] for n in nodeslist]
        y = [n.state[2] for n in nodeslist]
        scatter!(roomPlot, x,y, color=:black)

        edgeXs, edgeYs = [], []
        for e in edgeslist
            startN = e.startNode
            endN = e.endNode
            x1,y1 = startN.state[1], startN.state[2]
            x2,y2 = endN.state[1], endN.state[2]
            push!(edgeXs, x1, x2, NaN) #the NaNs, keep spaces between edges correctly unplotted
            push!(edgeYs, y1, y2, NaN)
        end

        plot!(roomPlot, edgeXs, edgeYs, color=:tan, linewidth=1)

        rrtPlot = plotSolPath(roomPlot, solPath)

        title!(rrtPlot, title)
        plot!(rrtPlot, legend=false, size=(600,600), xaxis=((-2,22), 0:1:20 ),
              yaxis=((-2,22), 0:1:20), foreground_color_grid= :black, titlefont=font("Arial", 10))

        return rrtPlot
    end

    function plotSolPath(aPlot, solPath)
        #print("\n ---Solution Path----- \n")
        if solPath != Void
            xPath = [n.state[1] for n in solPath]
            yPath = [n.state[2] for n in solPath]
            xstart, ystart = solPath[1].state

            xend, yend = solPath[end].state
            scatter!(aPlot, [xstart], [ystart],
                     markercolor= :red, markershape = :circle,  markersize = 6, markerstrokealpha = 0.5, markerstrokewidth=1)
            scatter!(aPlot, [xend], [yend],
                     markerstrokecolor = :green, markershape = :dtriangle,  markersize = 5, markerstrokealpha = 1, markerstrokewidth=5)
            plot!(aPlot, xPath, yPath, color = :orchid, linewidth=3, fillalpha = 0.3)
        else
           # #print("\n --- No solution path found ----- \n")
        end

        return aPlot
    end

end


text:
module CSV

const TypeVec = Vector{DataType}
# ============================================================================ #
struct Row{T<:AbstractString}
    data::T
    delim::Char
end
# ---------------------------------------------------------------------------- #
Base.iterate(r::Row) = iterate(r, 1)
# ---------------------------------------------------------------------------- #
function Base.iterate(r::Row, klast::Int)
    klast < 1 && return nothing

    while isspace(r.data[klast])
        klast += 1
    end

    if r.data[klast] == '"'
        # find the next '"' that is not escaped (preceded by a '\')
        # k = search(r.data, '"', klast+1)
        k = findnext(isequal('"'), r.data, klast+1)
        while r.data[k-1] == '\\'
            # k = search(r.data, '"', k+1)
            k = findnext(isequal('"'), r.data, k+1)
        end

        # k points to the char after the closing '"'
        k += 1

        # knext needs to point to the char after the closing r.delim or 0
        # knext = search(r.data, r.delim, k)
        knext = findnext(isequal(r.delim), r.data, k)
        knext = knext == nothing ? 0 : knext + 1
    else
        # element is not quoted, simply find the next delim
        # k = search(r.data, r.delim, klast)
        k = findnext(isequal(r.delim), r.data, klast)
        if k == nothing
            knext = 0
            k = length(r.data) + 1
        elseif k == length(r.data)
            # string ends with delim (i.e. last element of this row is empty)
            knext = 0
        else
            knext = k + 1
        end
    end
    return SubString(r.data, klast, k-1), knext
end
# ---------------------------------------------------------------------------- #
# Base.done(r::Row, k::Int) = k < 1
Base.IteratorSize(r::Row) = SizeUnknown()
Base.IteratorEltype(r::Row) = HasEltype()
# ---------------------------------------------------------------------------- #
Base.eltype(r::Row{T}) where T<:AbstractString = T
# ---------------------------------------------------------------------------- #
function Base.collect(r::Row)
    out = Vector{String}()
    for x in r
        push!(out, x)
    end
    return out
end
# ============================================================================ #
function parse(ifile::String, delim::Char=',', types::TypeVec=TypeVec())
    return try
        if isfile(ifile)
            return parse_impl(read(ifile, String), delim, types)
        elseif isdir(dirname(ifile))
            error("Input appears to be an invalid file path...")
        else
            return parse_impl(ifile, delim, types)
        end
    catch err
        if isa(err, Base.IOError)
            return parse_impl(ifile, delim, types)
        else
            rethrow(err)
        end
    end
end
# ---------------------------------------------------------------------------- #
function parse_impl(csvstr::String, delim::Char=',', types::TypeVec=TypeVec())
    lines = split(csvstr, r"\n|\r\n|\r", limit=0, keepempty=false)
    hdr = split_line(lines[1], delim)
    if isempty(types)
        types = fill(Any, length(hdr))
    elseif length(types) < length(hdr)
        append!(types, fill(Any, length(hdr) - length(types)))
    end

    data = Dict{String, Vector}(
        (hdr[k] => Vector{types[k]}(undef, length(lines)-1) for k in 1:length(hdr))
    )

    for k = 2:length(lines)
        j = 1
        for item in Row(lines[k], delim)
            data[hdr[j]][k-1] = parse_item(types[j], item)
            j += 1
        end
    end
    return data
end
# ============================================================================ #
isquoted(str::AbstractString) = occursin(r"^\"((?:[^\"\\]|\\.)*)\"$", strip(str))
# ============================================================================ #
dequote(str::AbstractString) = replace(str, r"(?:^\s*\")|(?:\"\s*$)"=>"")
# ============================================================================ #
function split_line(str::AbstractString, delim::Char)
    out = Vector{String}()
    for x in Row(str, delim)
        push!(out, dequote(x))
    end
    return out
end
# ============================================================================ #
parse_item(::Type{T}, item::AbstractString) where {T<:Number} =
    Base.parse(T, strip(item))
# ============================================================================ #
function parse_item(::Type, item::AbstractString)
    return isquoted(item) ? dequote(item) : strip(item)
end
# ============================================================================ #
end


text:
module TestSoilGrid

using Base.Test
using SoilTracers.SoilGrid

@time @testset "SoilGrid constructor" begin
    required_fields =
        Set(Symbol[:grid_size, :grid_top, :grid_bottom, :grid_node])

    # test case: empty soil grid
    grid1 = FVGrid()
    @test grid1.level > 0
    @test grid1.top < grid1.bottom
    @test isa(grid1.texture, String)
    @test Set(grid1.profiles.colindex.names) >= required_fields
    @test all(grid1.profiles[:grid_top] .< grid1.profiles[:grid_node])
    @test all(grid1.profiles[:grid_node] .< grid1.profiles[:grid_bottom])

    # test case: a soil grid with more fields
    grid2 = FVGrid(
        25, 0., 1., "sandy loam", [],
        ["temp", "moisture", "porosity", "pH", "orgC"],
        ["co2"])
    @test grid2.level > 0
    @test grid2.top < grid2.bottom
    @test isa(grid2.texture, String)
    @test Set(grid2.profiles.colindex.names) >= required_fields
    @test all(grid2.profiles[:grid_top] .< grid2.profiles[:grid_node])
    @test all(grid2.profiles[:grid_node] .< grid2.profiles[:grid_bottom])
    @test all([
        in(field, grid2.profiles.colindex.names) for
        field in Symbol.(["temp", "moisture", "porosity", "pH", "orgC"])])

    # test case: ensure that top depth < bottom depth always true
    grid3 = FVGrid(20, 1.5, 0.)
    @test grid3.level > 0
    @test grid3.top < grid3.bottom
end

@time @testset "SoilGrid methods" begin
    grid = FVGrid(
        25, 0., 1., "sandy loam", [],
        ["temp", "moisture", "porosity", "pH", "orgC"],
        ["co2"])

    # getprofile
    co2_profile = getprofile(grid, :co2)
    @test typeof(co2_profile) <: AbstractArray
    @test length(co2_profile) == grid.level
    @test_throws ArgumentError getprofile(grid, :cos)

    # setprofile!
    co2 = 400.  # arbitrary unit
    setprofile!(grid, :co2, co2)
    co2_profile = getprofile(grid, :co2)
    @test all(isapprox.(co2_profile, co2))
    @test_throws ArgumentError setprofile!(grid, :cos, co2)

    # addprofile!
    cos = 500.  # arbitrary unit
    @test_throws ArgumentError addprofile!(grid, :co2)
    addprofile!(grid, :cos)
    @test :cos in grid.profiles.colindex.names

    # delprofile!
    delprofile!(grid, :cos)
    @test !(:cos in grid.profiles.colindex.names)
    @test_throws ArgumentError delprofile!(grid, :cos)

    # validategrid
    @test validategrid(grid) === nothing
    delprofile!(grid, :temp)
    @test_throws ArgumentError validategrid(grid)
end

end


text:
using SpecialFunction

clipped_intensity(edge_pos, beam_radius) = 1 / 2 * (1 - erf(sqrt(2) * edge_pos / beam_radius))

export clipped_intensity


text:
macro staticvar(init)
    var = gensym()
    __module__.eval(:(const $var = $init))
    var = esc(var)
    quote
        global $var
        $var
    end
end

macro staticdef(ex)
    @capture(ex, name_::T_ = val_) || error("invalid @staticvar")
    ref = Ref{__module__.eval(T)}()
    set = Ref(false)
    :($(esc(name)) = if $set[]
        $ref[]
    else
        $ref[] = $(esc(ex))
        $set[] = true
        $ref[]
    end)
end

unix2date(t) = Date(unix2datetime(t))
unix2time(t) = Time(unix2datetime(t))
unix2hour(x) = x % (24 * 3600) / 3600

unix2str8(t) = Dates.format(unix2datetime(t), dateformat"yyyymmdd")
unix2str6(t) = Dates.format(unix2datetime(t), dateformat"yymmdd")
str2date(str) = Date(replace(str, '-' => ""), dateformat"yyyymmdd")
str2datetime(str) = DateTime(replace(str, '-' => ""), dateformat"yyyymmdd")
str2unix(str) = datetime2unix(str2datetime(str))

unix2int(t) = parse(Int, unix2str8(t))
int2unix(i) = str2unix(string(i))

function mmap_tempname()
    isdir(".mempool") || mkdir(".mempool")
    file = joinpath(".mempool", randstring())
end

function mcopy(x::AbstractArray{T, N}) where {T, N}
    src = mmap_tempname()
    fid = open(src, "w+")
    write(fid, x)
    seekstart(fid)
    xm = Mmap.mmap(fid, Array{T, N}, size(x))
    close(fid)
    finalizer(z -> rm(src), xm)
    return xm
end

mcopy(x::Series) = Series(mcopy(x.values), index = x.index, name = x.name)

mcopy(x::DataFrame) = DataFrame(mcopy(x.values), index = x.index, columns = x.columns)

part(df::DataFrame) = nprocs() > 1 ? part(df.iloc, 1) : df

lngstconsec(s) = (!s).cumsum()[s].value_counts().max()

function reindex_columns(df, columns)
    itrs = [df[c].sort_values().unique() for c in columns]
    multidx = pd.MultiIndex.from_product(itrs, names = columns)
    df.set_index(columns, inplace = true)
    df = df.reindex(index = multidx)
    df.reset_index(inplace = true)
    return df
end

function concat_hdfs(dst, srcs; key = "df")
    store = pd.HDFStore(dst, "w")
    try
        for src in srcs
            df = pd.read_hdf(src, key)
            store.append(key, df, index = false)
        end
    finally
        store.close()
    end
    return dst
end

hasheader(src) = length(collect(eachmatch(r"\d+\.\d+", readline(src)))) <= 2

function concat_txts(dst, srcs)
    header = hasheader(srcs[1])
    open(dst, "w") do fid
        header && println(fid, readline(srcs[1]))
        @showprogress 10 dst for src in srcs
            open(src, "r") do f
                header && readline(f)
                write(fid, read(f))
            end
        end
    end
    return dst
end

function concat_lagged(df::DataFrame, by; maxlag = 10)
    df = pd.concat([df.groupby(by).shift(l).dropna() for l in 0:(maxlag - 1)], axis = 1)
    df.columns = [string(c, '_', l) for l in 0:(maxlag - 1) for c in df.columns]
    df.reset_index(inplace = true, drop = true)
    return df
end

function concat_lagged(x; maxlag = 10)
    F, N, T = size(x)
    x′ = fill!(similar(x, maxlag * F, N, T), 0)
    for t in 1:T, n in 1:N, f in 1:F
        s = x[f, n, t]
        for t′ in t:min(t + maxlag - 1, T)
            x′[f + (t′ - t) * F, n, t′] = s
        end
    end
    return x′
end

function read_hdf5(src; columns = nothing, mmaparrays = true)
    df = DataFrame()
    h5open(src, "r+") do fid
        columns = something(columns, names(fid))
        @showprogress 10 "read_hdf5..." for c in columns ∩ names(fid)
            dset = fid[c]
            if Sys.iswindows() || !mmaparrays || !ismmappable(dset)
                df[c] = read(dset)
            else
                df[c] = readmmap(dset)
            end
        end
    end
    return df
end

function to_hdf5(df, dst)
    isfile(dst) && rm(dst)
    h5open(dst, "w") do fid
        @showprogress 10 "to_hdf5..." for c in df.columns
            fid[c] = Array(df[c])
        end
    end
    return dst
end

to_dict(x) = Dict{Symbol, Any}(s => getfield(x, s) for s in fieldnames(typeof(x)))

to_struct(::Type{T}, d) where T = T([d[s] for s in fieldnames(T)]...)

divavg!(x) = (x ./= mean(x))

function rollindices(ti, tf, Δtb, Δtf)
    [(string(t - Δtb):string(t - Day(1)),
    string(t):string(t + Δtf - Day(1)))
    for t in (ti + Δtf):Δtf:tf]
end

nunique(x) = length(Set(x))

function sortednunique(f, x)
    n = 1
    yl = f(x[1])
    @inbounds for xi in x
        yi = f(xi)
        n = ifelse(yi != yl, n + 1, n)
        yl = yi
    end
    return n
end

sortednunique(x) = sortednunique(identity, x)

function sortedunique(f, x)
    n = 1
    yl = f(x[1])
    yu = [yl]
    @inbounds for xi in x
        yi = f(xi)
        if yi != yl
            n += 1
            push!(yu, yi)
        end
        yl = yi
    end
    return yu
end

sortedunique(x) = sortedunique(identity, x)

Mmap.sync!(m::AbstractArray) = isdefined(m, :parent) && Mmap.sync!(m.parent)

parsefreq(freq) = freq

function parsefreq(freq::String)
    @eval let S = 1, T = 60, H = 60T, D = 24H, M = 30D, Y = 12M
        $(Meta.parse(freq))
    end
end

idxmap(x) = Dict(zip(x, axes(x, 1)))

⧶(x, y) = x / y
⧶(x, y::AbstractFloat) = x / (y + eps(y))
⧶(x, y::Integer) = ifelse(x == y == 0, zero(x), x / y)

function from_category(sr, T = MLString{8})
    sr = sr.astype("category")
    cats = Array(sr.cat.categories)
    codes = (sr.cat.codes + 1).values
    T[cats[i] for i in codes]
end

function to_category(x)
    sr = pd.Series(x).astype("category")
    sr.cat.categories = sr.cat.categories.astype("str")
    return sr
end

if PandasLite.version() >= v"0.25"
    function to_category(x::AbstractArray{<:MLString})
        sr = pd.Series(x).astype("category")
        sr.cat.categories = sr.cat.categories.str.decode("utf-8", "ignore")
        return sr
    end
end

Base.dropdims(f, A::AbstractArray; dims) = dropdims(f(A, dims = dims), dims = dims)

function arr2rng(x)
    isempty(x) && return x
    r = UnitRange(extrema(x)...)
    r == x ? r : x
end

isna(x) = iszero(x) | isnan(x)

fillnan(x) = ifelse(isnan(x), zero(x), x)

fillnan!(x) = map!(fillnan, x, x)


text:
"""
    correlationmatrix(hopmatrix, n_occupied)

calculates the two-body correlation matrix lambda ``Λ = ⟨a^†_i a_j⟩``
by full diagonalization of given Hamiltonian matrix.

Here, the correlation matrix Λ is made from its diagonal form Γ, which
has n_occupied 1s--corresponding to the lowest eigenvalues--and the
rest of the diagonal are 0s. We have ``Λ = U^{*} Γ U^{T}``

"""
function correlationmatrix(hopmatrix::Matrix{Float64},
                           n_occupied ::Int64)

    lx = size(hopmatrix)[1]
    @assert 0 < n_occupied && n_occupied <= lx

    fact = eigen(Symmetric(hopmatrix))
    println(fact.values[n_occupied-2:n_occupied+2])
    e_fermi = fact.values[n_occupied]
    n=1
    while fact.values[n_occupied+n] - e_fermi < 1.e-10
        n=n+1
    end
    n > 1 && @warn "Hamiltonian is $n-fold degenerate at the given filling factor $n_occupied/$lx"
    vecs = fact.vectors[:,1:n_occupied]
    #vecs = eigvecs(Symmetric(hopmatrix))[:, 1:n_occupied]
    return vecs * transpose(vecs)
end

function correlationmatrix(hopmatrix::Matrix{ComplexF64},
                           n_occupied ::Int64)

    lx = size(hopmatrix)[1]
    @assert 0 < n_occupied && n_occupied <= lx

    fact = eigen(Hermitian(hopmatrix))
    e_fermi = fact.values[n_occupied]
    n=1
    while fact.values[n_occupied+n] - e_fermi < 1.e-10
        n=n+1
    end
    n > 1 && @warn "Hamiltonian is $n-fold degenerate at the given filling factor $n_occupied/$lx"
    vecs = fact.vectors[:,1:n_occupied]
    #vecs = eigvecs(Hermitian(hopmatrix))[:, 1:n_occupied]
    return conj(vecs) * transpose(vecs)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
# ---
# title: 318. Maximum Product of Word Lengths
# id: problem318
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/maximum-product-of-word-lengths/description/>
# hidden: true
# ---
#
# Given a string array `words`, find the maximum value of `length(word[i]) *
# length(word[j])` where the two words do not share common letters. You may
# assume that each word will contain only lower case letters. If no such two
# words exist, return 0.
#
# **Example 1:**
#
#
#
#     Input: ["abcw","baz","foo","bar","xtfn","abcdef"]
#     Output:16
#     Explanation:The two words can be "abcw", "xtfn".
#
# **Example 2:**
#
#
#
#     Input: ["a","ab","abc","d","cd","bcd","abcd"]
#     Output:4
#     Explanation:The two words can be "ab", "cd".
#
# **Example 3:**
#
#
#
#     Input: ["a","aa","aaa","aaaa"]
#     Output:0
#     Explanation:No such pair of words.
#
#
#
#
# **Constraints:**
#
#   * `0 <= words.length <= 10^3`
#   * `0 <= words[i].length <= 10^3`
#   * `words[i]` consists only of lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
@testset "Scalar spherical harmonic mode indices" begin
    lmax = 10
    N = lmax + 1
    M = 2 * N - 1

    mode_seen = zeros(Int, N, M)
    for l in 0:lmax, m in (-l):l
        mode_seen[sph_mode(l, m)] += 1
    end
    @test mode_seen == sphones(Bool, N, M)
end

@testset "Scalar spherical harmonics: simple modes ($T)" for T in [Float64,
                                                              Complex{Float64}]
    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    # Test l=0, m=0
    F = T[Y_0_0(θ, ϕ) for θ in Θ, ϕ in Φ]
    C = sph_transform(F)
    @test C ≈ unit(sph_mode(0, 0), size(C))
    F′ = sph_evaluate(C)
    @test F′ ≈ F

    # Test l=1, m=0
    F = T[Y_1_0(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test size(F) == (N, M)
    C = sph_transform(F)
    @test C ≈ unit(sph_mode(1, 0), size(C))
    F′ = sph_evaluate(C)
    @test F′ ≈ F

    # Test l=1, m=1
    F = T[Y_1p1(θ, ϕ) for θ in Θ, ϕ in Φ]
    C = sph_transform(F)
    @test C ≈ unit(sph_mode(1, 1), size(C))
    F′ = sph_evaluate(C)
    @test F′ ≈ F

    # Test l=1, m=-1
    F = T[Y_1m1(θ, ϕ) for θ in Θ, ϕ in Φ]
    C = sph_transform(F)
    @test C ≈ unit(sph_mode(1, -1), size(C))
    F′ = sph_evaluate(C)
    @test F′ ≈ F
end

@testset "Scalar spherical harmonics: simple modes ($T)" for T in [Float64,
                                                              Complex{Float64}]
    lmax = 10

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    lmax_test = 4
    for l in 0:lmax_test, m in (-l):l
        F = T[Ylm(l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        C = sph_transform(F)
        @test C ≈ unit(sph_mode(l, m), size(C))
        F′ = sph_evaluate(C)
        @test F′ ≈ F
    end
end

@testset "Scalar spherical harmonics: linearity ($T)" for T in [Float64,
                                                           Complex{Float64}]
    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    F = zeros(T, N, M)
    C = sph_transform(F)
    @test iszero(C)

    F1 = randn(T, N, M)
    F2 = randn(T, N, M)
    C1 = sph_transform(F1)
    C2 = sph_transform(F2)
    F = F1 + F2
    C = sph_transform(F)
    @test C ≈ C1 + C2

    F = randn(T, N, M)
    α = randn(T)
    C = sph_transform(F)
    Cα = sph_transform(α * F)
    @test Cα ≈ α * C
end

@testset "Scalar spherical harmonics: duality ($T)" for T in [Float64,
                                                         Complex{Float64}]
    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    F = randn(T, N, M)
    C = sph_transform(F)
    @test size(C) == size(F)
    @test eltype(C) == eltype(F)
    F′ = sph_evaluate(C)
    @test F′ ≈ F

    C = randn(T, N, M)
    F = sph_evaluate(C)
    @test size(F) == size(C)
    @test eltype(F) == eltype(C)
    C′ = sph_transform(F)
    @test C′ ≈ C
end

@testset "Scalar spherical harmonics: orthonormality ($T)" for T in [Float64,
                                                                Complex{Float64}]
    lmax = 100
    atol = 4 / lmax^2

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    lmax_test = 4
    for l in 0:lmax_test, m in (-l):l
        C = zeros(T, N, M)
        C[sph_mode(l, m)] = 1
        F = sph_evaluate(C)
        for l′ in 0:lmax_test, m′ in (-l′):l′
            C′ = zeros(T, N, M)
            C′[sph_mode(l′, m′)] = 1
            F′ = sph_evaluate(C′)
            int = integrate(F, F′)
            δ = l == l′ && m == m′
            @test isapprox(int, δ; atol=atol)
        end
    end
end

@testset "Scalar spherical harmonics: Laplacian ($T)" for T in [Float64,
                                                           Complex{Float64}]
    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    # z + 2x
    F = T[cos(θ) + 2 * sin(θ) * cos(ϕ) for θ in Θ, ϕ in Φ]
    C = sph_transform(F)
    ΔC = sph_laplace(C)
    ΔF = sph_evaluate(ΔC)
    @test ΔF ≈ -2F
end


text:
function exchangeEffectiveField(coupling::Matrix, u)
    # coupling is a NxN array which describes which sites couple to which
    # returns an array with effective field
    Heff = zeros(Float64, size(coupling)[1], 3)

    for i in 1:size(coupling)[1]
        for j in 1:size(coupling)[1]
            @. Heff[i,:] += coupling[i,j] .* u[j,:]
        end
    end
    Heff
end

function exchangeEffectiveField!(Heff, coupling::Matrix, u)
    # coupling is a NxN array which describes which sites couple to which
    # returns an array with effective field
    for i in 1:size(coupling)[1]
        for j in 1:size(coupling)[1]
            uj = @view u[j,:]
            cij = @view coupling[i,j]
            @. Heff[i,:] += cij .* uj
        end
    end
end

function DMIEffectiveField(coupling::Matrix, Dfield::Vector, u)
    # coupling is a NxN array which describes which sites couple to which
    # returns an array with effective field
    Heff = zeros(Float64, size(coupling)[1], 3)

    for i in 1:size(coupling)[1]
        for j in 1:size(coupling)[1]
            Heff[i,:] .-= coupling[i,j] * cross(Dfield, u[j,:])
        end
    end
    Heff
end

function magneticEffectiveField(Hfield::Vector, N)
    Heff = zeros(Float64, N, 3)
    for i in 1:N
        @. Heff[i,:] += Hfield
    end
    return Heff
end

function magneticEffectiveField!(Heff, Hfield::Vector, N)
    for i in 1:N
        @. Heff[i,:] += Hfield
    end
end

function spinTorqueEffectiveField(Js, u)
    Heff = zeros(Float64, size(Js)[1], 3)
    for i in 1:size(Js)[1]
        Heff[i,:] = cross(Js[i,:], u[i,:])
    end
    return Heff
end

# function generateEOM(Heff, α, vars::VarsCart)
#     res = zeros(Num, length(vars.x), 3)
#     for i in 1:length(vars.x)
#         S = [vars.x[i], vars.y[i], vars.z[i]]
#         res[i, :] = - cross(S, Heff[i,:]) - α * cross(S, cross(S, Heff[i,:]))
#     end
#     D = Differential(vars.t)
#     return [D.(vars.x) .~ res[:, 1]; D.(vars.y) .~ res[:, 2]; D.(vars.z) .~ res[:, 3]]
# end

function initRandomSpin(N)
    # boltzmann is proportional to θ (for now)
    # the angle ϕ is chosen randomly between 0 and 2pi
    ϕ = rand(N) * 2π
    θ = rand(N) * π
    mx = sin.(θ) .* cos.(ϕ)
    my = sin.(θ) .* sin.(ϕ)
    mz = cos.(θ)
    m0 = zeros(N, 3)
    m0[:,1] = mx
    m0[:,2] = my
    m0[:,3] = mz
    return m0
end

function initBoltzmann(N, λ)
    # boltzmann is proportional to θ (for now)
    # the angle ϕ is chosen randomly between 0 and 2pi
    dist = Exponential(λ / 2)
    θ = rand(dist, N)
    ϕ = rand(Float64, N) * 2π
    mx = sin.(θ) .* cos.(ϕ)
    my = sin.(θ) .* sin.(ϕ)
    mz = cos.(θ)
    m0 = zeros(N, 3)
    m0[:,1] = mx
    m0[:,2] = my
    m0[:,3] = mz
    return m0

end

max_stars_repo_name: pgunnink/LLGTools

max_issues_repo_name: pgunnink/LLGTools

max_forks_repo_name: pgunnink/LLGTools


text:
using PLSRegressor
using Test

include("./pls1_test.jl")
include("./pls2_test.jl")
include("./utils_test.jl")
include("./kpls_test.jl")


text:
using MAT
include("../src/hypergraphs.jl")
include("../src/SparseCard.jl")

## Load

datasets = ["020_smallplant","001_oct"]

for numsuperpix = [500; 200]

    if numsuperpix == 200
        lambdas = [0.3; 0.1]
        lambda2 = 0.002

    elseif numsuperpix == 500
        lambdas = [0.08; 0.01]

        lambda2 = 0.008
    end

for jj = 1:2

dataset = datasets[jj]
lambda1 = lambdas[jj]
M = matread("../data/data_reflections/$(dataset)_p$numsuperpix.mat")

# The graph is a grid: each node is a pixel in a picture
# This gives superpixel regions
sup_img = M["sup_img"]

# Load
M2 = matread("../data/data_reflections/$(dataset)_wts.mat")
W1 = M2["W1"] # edge weights for vertical edges
W2 = M2["W2"] # edge weights for horizontal edges
unary = M2["unary"] # this is the s-t edge weights, already linearized

# number of columns
col = size(sup_img,2)

# number of rows
row = size(sup_img,1)

# go from (i,j) node ID to linear node ID
node = (i,j)->i + (j-1)*row

## Build edge information
I = Vector{Int64}()
J = Vector{Int64}()
W = Vector{Float64}()

# Type 1 edges: terminal edges
tvec = abs.((unary .> 0).* unary)
svec = abs.((unary .< 0).* unary)

# Type 2 edges: vertical edges, organized into rows
for i = 1:row-1
    for j = 1:col
        v1 = node(i,j)      # first node in edge
        v2 = node(i+1,j)    # second node in edge
        push!(I,v1)
        push!(J,v2)

        # Use edge weights as defined by Jegelka et al.
        push!(W,lambda1*W1[i,j])
    end
end

# Type 3 edges: horizontal edges, organized into columns
for i = 1:col-1
    for j = 1:row
        v1 = node(j,i)      # first node in edge
        v2 = node(j,i+1)    # second node in edge
        push!(I,v1)
        push!(J,v2)

        # Use edge weights as defined by Jegelka et al.
        push!(W,lambda1*W2[j,i])
    end
end
N = length(unary)
A = sparse(I,J,W,N,N)
A = A+sparse(A')
matwrite("Graph_noclique_$(dataset)_$numsuperpix.mat",Dict("A"=>A,"svec"=>svec, "tvec"=>tvec))

end
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM