text:
using JuBrain
using Base.Test

model=NetworkModel()

model.parameters="""
τᵥ = 20*ms;  τₑ = 5*ms;   τᵢ = 10*ms
Vt = -50*mV; Vᵣ = -60*mV; El = -49*mV"""
model.equations="dv/dt=(gₑ+gᵢ-(v-El))/τᵥ;dgₑ/dt=-gₑ/τₑ;dgᵢ/dt=-gᵢ/τᵢ"
model.spike="v.>Vt"
model.reset="v=Vᵣ"
model.synapse="gₑ+=wₑ;gᵢ+=wᵢ"""

model.groups=makeGroups(groupName=["e","i"],groupSize=[3200,800])
initialize(model,groups="all",expr="v=Vᵣ+rand()*(Vt-Vᵣ)")
spConnect(model,pre="e",post="all",expr="wₑ=(60*0.27/10)*mV",p=0.02)
spConnect(model,pre="i",post="all",expr="wᵢ=(-20*4.5/10)*mV",p=0.02)

buildNetwork(name="temp",model=model,duration=1000)


text:
# TODO: Adapt python model below to new Julia format
# ################################################################################
# ##                          crystalline MgF2                                  ##
# ################################################################################
# def n_MgF2_sym(axis='e'):
#     """This function creates a symbolic representation (using SymPy) of the
#     Sellmeier Equation model for the wavelength dependence
#     of crystalline MgF2's ordinary and extraordinary indices of refraction as
#     well as the index of thin-film amorphous MgF2.
#     Sellmeier coefficients for crystalline MgF2 are taken from
#     "Refractive properties of magnesium fluoride"
#     by Dodge, Applied Optics 23 (12), pp.1980-1985 (1984)
#     https://doi.org/10.1364/AO.23.001980
#     Sellmeier coefficients for amorphous MgF2 are taken from
#     "Self-consistent optical constants of MgF2, LaF3, and CeF3 films"
#     by Rodríguez-de Marcos, et al. Optical Materials Express 7 (3) (2017)
#     https://doi.org/10.1364/OME.7.000989
#     This model is then exported to other functions that use it and its
#     derivatives to return index, group index and GVD values as a function
#     of temperature and wavelength.
#     Variable units are lm in [um] and T in [deg C]
#     """
#     if axis is 'e':
#         coeffs = [  (0.41344023,0.03684262), # A_i,λ_i Lorentz oscillator strength [1] and resonance wavelength [μm]
#                     (0.50497499,0.09076162),
#                     (2.4904862,23.771995),
#                     ]
#     elif axis is 'o':
#         coeffs = [  (0.48755108,0.04338408), # A_i,λ_i Lorentz oscillator strength [1] and resonance wavelength [μm]
#                     (0.39875031,0.09461442),
#                     (2.3120353,23.793604),
#                     ]
#     elif axis is 'a':
#         coeffs = [  (1.73,.0805), # A_i,λ_i Lorentz oscillator strength [1] and resonance wavelength [μm]
#                     ]
#     else:
#         raise Exception('unrecognized axis! must be "e", "o" or "a"')
#     lm, T = sp.symbols('lm T')
#     A0,λ0 = coeffs[0]
#     oscillators = A0 * lm**2 / (lm**2 - λ0**2)
#     if len(coeffs)>1:
#         for Ai,λi in coeffs[1:]:
#             oscillators += Ai * lm**2 / (lm**2 - λi**2)
#     n_sym = sp.sqrt( 1 + oscillators )
#     return lm, T, n_sym
#
# def n_MgF2(lm_in,T_in=300*u.degK,axis='e'):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the index of refraction of crystalline (axis='o' or 'e') and
#     amorphous (axis='a') MgF2.
#     Sellmeier coefficients for crystalline MgF2 are taken from
#     "Refractive properties of magnesium fluoride"
#     by Dodge, Applied Optics 23 (12), pp.1980-1985 (1984)
#     https://doi.org/10.1364/AO.23.001980
#     Sellmeier coefficients for amorphous MgF2 are taken from
#     "Self-consistent optical constants of MgF2, LaF3, and CeF3 films"
#     by Rodríguez-de Marcos, et al. Optical Materials Express 7 (3) (2017)
#     https://doi.org/10.1364/OME.7.000989
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = np.array([Q_(lm_in).to(u.um).magnitude]).flatten()
#     T_C = np.array([Q_(T_in).to(u.degC).magnitude]).flatten()
#     lm, T, n_sym = n_MgF2_sym(axis=axis)
#     n = sp.lambdify([lm,T],n_sym,'numpy')
#     output = np.zeros((T_C.size,lm_um.size))
#     for T_idx, TT in enumerate(T_C):
#         output[T_idx,:] = n(lm_um, T_C[T_idx])
#     return output.squeeze()
#
# def n_g_MgF2(lm_in,T_in=300*u.degK,axis='e'):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group index of refraction of crystalline (axis='o' or 'e') and
#     amorphous (axis='a') MgF2.
#     Sellmeier coefficients for crystalline MgF2 are taken from
#     "Refractive properties of magnesium fluoride"
#     by Dodge, Applied Optics 23 (12), pp.1980-1985 (1984)
#     https://doi.org/10.1364/AO.23.001980
#     Sellmeier coefficients for amorphous MgF2 are taken from
#     "Self-consistent optical constants of MgF2, LaF3, and CeF3 films"
#     by Rodríguez-de Marcos, et al. Optical Materials Express 7 (3) (2017)
#     https://doi.org/10.1364/OME.7.000989
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_MgF2_sym(axis=axis)
#     n_sym_prime = sp.diff(n_sym,lm)
#     n_g_sym = n_sym - lm*n_sym_prime
#     n_g = sp.lambdify([lm,T],n_g_sym,'numpy')
#     return n_g(lm_um, T_C)
#
# def gvd_MgF2(lm_in,T_in=300*u.degK,axis='e'):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group velocity dispersion of crystalline (axis='o' or 'e') and
#     amorphous (axis='a') MgF2.
#     Sellmeier coefficients for crystalline MgF2 are taken from
#     "Refractive properties of magnesium fluoride"
#     by Dodge, Applied Optics 23 (12), pp.1980-1985 (1984)
#     https://doi.org/10.1364/AO.23.001980
#     Sellmeier coefficients for amorphous MgF2 are taken from
#     "Self-consistent optical constants of MgF2, LaF3, and CeF3 films"
#     by Rodríguez-de Marcos, et al. Optical Materials Express 7 (3) (2017)
#     https://doi.org/10.1364/OME.7.000989
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_MgF2_sym(axis=axis)
#     n_sym_double_prime = sp.diff(n_sym,lm,lm)
#     c = Q_(3e8,'m/s') # unitful definition of speed of light
#     gvd_sym_no_prefactors = (lm**3)*n_sym_double_prime # symbolic gvd without unitful prefactors, to be made unitful below
#     gvd_no_prefactors = sp.lambdify([lm,T],gvd_sym_no_prefactors,'numpy') # numerical gvd without unitful prefactors
#     gvd = (1 / (2 * np.pi * (c**2))) * Q_(gvd_no_prefactors(lm_um,T_C),'um')
#     return gvd.to('fs**2 / mm')
#
#
#

ext: jl


text:
const ASICAMERA_ID_MAX = 128

# Type and error enumerations

@cenum ASI_CONTROL_TYPE::UInt32 begin
    ASI_GAIN = 0
    ASI_EXPOSURE = 1  # Exposure time in μs
    ASI_GAMMA = 2     # 1-100, default: 50
    ASI_WB_R = 3      # white balance, red component
    ASI_WB_B = 4      # white balance, blue
    ASI_OFFSET = 5    # pixel value offset / bias
    ASI_BANDWIDTHOVERLOAD = 6  # data transfer rate percentage
    ASI_OVERCLOCK = 7
    ASI_TEMPERATURE = 8  # 10 times the actual temperature
    ASI_FLIP = 9
    ASI_AUTO_MAX_GAIN = 10  # Maximum gain when auto adjust
    ASI_AUTO_MAX_EXP = 11   # Maximum exposure time when auto adjust, μs
    ASI_AUTO_TARGET_BRIGHTNESS = 12  # Target brightness when auto adjust
    ASI_HARDWARE_BIN = 13
    ASI_HIGH_SPEED_MODE = 14
    ASI_COOLER_POWER_PERC = 15
    ASI_TARGET_TEMP = 16  # °C, don't multiply by 10
    ASI_COOLER_ON = 17
    ASI_MONO_BIN = 18  # lead to a smaller grid at software bin mode for color camera?!
    ASI_FAN_ON = 19
    ASI_PATTERN_ADJUST = 20  # currently only supported by 1600 mono camera
    ASI_ANTI_DEW_HEATER = 21
end

const ASI_BRIGHTNESS = ASI_OFFSET
const ASI_AUTO_MAX_BRIGHTNESS = ASI_AUTO_TARGET_BRIGHTNESS

@cenum ASI_BAYER_PATTERN::UInt32 begin
    ASI_BAYER_RG = 0
    ASI_BAYER_BG = 1
    ASI_BAYER_GR = 2
    ASI_BAYER_GB = 3
end

@cenum ASI_IMG_TYPE::Int32 begin
    ASI_IMG_RAW8 = 0
    ASI_IMG_RGB24 = 1
    ASI_IMG_RAW16 = 2
    ASI_IMG_Y8 = 3
    ASI_IMG_END = -1
end

@cenum ASI_GUIDE_DIRECTION::UInt32 begin
    ASI_GUIDE_NORTH = 0
    ASI_GUIDE_SOUTH = 1
    ASI_GUIDE_EAST = 2
    ASI_GUIDE_WEST = 3
end

@cenum ASI_FLIP_STATUS::UInt32 begin
    ASI_FLIP_NONE = 0
    ASI_FLIP_HORIZ = 1
    ASI_FLIP_VERT = 2
    ASI_FLIP_BOTH = 3
end

@cenum ASI_CAMERA_MODE::Int32 begin
    ASI_MODE_NORMAL = 0
    ASI_MODE_TRIG_SOFT_EDGE = 1
    ASI_MODE_TRIG_RISE_EDGE = 2
    ASI_MODE_TRIG_FALL_EDGE = 3
    ASI_MODE_TRIG_SOFT_LEVEL = 4
    ASI_MODE_TRIG_HIGH_LEVEL = 5
    ASI_MODE_TRIG_LOW_LEVEL = 6
    ASI_MODE_END = -1
end

@cenum ASI_TRIG_OUTPUT::Int32 begin
    ASI_TRIG_OUTPUT_PINA = 0
    ASI_TRIG_OUTPUT_PINB = 1
    ASI_TRIG_OUTPUT_NONE = -1
end

const ASI_TRIG_OUTPUT_PIN = ASI_TRIG_OUTPUT

@cenum ASI_ERROR_CODE::UInt32 begin
    ASI_SUCCESS = 0
    ASI_ERROR_INVALID_INDEX = 1
    ASI_ERROR_INVALID_ID = 2
    ASI_ERROR_INVALID_CONTROL_TYPE = 3
    ASI_ERROR_CAMERA_CLOSED = 4
    ASI_ERROR_CAMERA_REMOVED = 5
    ASI_ERROR_INVALID_PATH = 6
    ASI_ERROR_INVALID_FILEFORMAT = 7
    ASI_ERROR_INVALID_SIZE = 8
    ASI_ERROR_INVALID_IMGTYPE = 9
    ASI_ERROR_OUTOF_BOUNDARY = 10
    ASI_ERROR_TIMEOUT = 11
    ASI_ERROR_INVALID_SEQUENCE = 12
    ASI_ERROR_BUFFER_TOO_SMALL = 13
    ASI_ERROR_VIDEO_MODE_ACTIVE = 14
    ASI_ERROR_EXPOSURE_IN_PROGRESS = 15
    ASI_ERROR_GENERAL_ERROR = 16
    ASI_ERROR_INVALID_MODE = 17
    ASI_ERROR_END = 18
end

struct ASIError <: Exception
    code::ASI_ERROR_CODE
end

Base.showerror(io::IO, err::ASIError) = print(io, err.code)

struct ASIWrapperError <: Exception
    message::String
end

Base.showerror(io::IO, err::ASIWrapperError) = print(io, err.message)

@cenum ASI_BOOL::UInt32 begin
    ASI_FALSE = 0
    ASI_TRUE = 1
end

ASI_BOOL(b::Bool) = b ? ASI_TRUE : ASI_FALSE

@cenum ASI_EXPOSURE_STATUS::UInt32 begin
    ASI_EXP_IDLE = 0
    ASI_EXP_WORKING = 1
    ASI_EXP_SUCCESS = 2
    ASI_EXP_FAILED = 3
end

# Structs
# kept mutable for now, has to be checked
# added constructors

mutable struct _ASI_CAMERA_INFO
    Name::NTuple{64, Cchar}
    CameraID::Cint
    MaxHeight::Clong
    MaxWidth::Clong
    IsColorCam::ASI_BOOL
    BayerPattern::ASI_BAYER_PATTERN
    SupportedBins::NTuple{16, Cint}
    SupportedVideoFormat::NTuple{8, ASI_IMG_TYPE}
    PixelSize::Cdouble
    MechanicalShutter::ASI_BOOL
    ST4Port::ASI_BOOL
    IsCoolerCam::ASI_BOOL
    IsUSB3Host::ASI_BOOL
    IsUSB3Camera::ASI_BOOL
    ElecPerADU::Cfloat
    BitDepth::Cint
    IsTriggerCam::ASI_BOOL
    Unused::NTuple{16, UInt8}
end

_ASI_CAMERA_INFO() = _ASI_CAMERA_INFO(
    NTuple{64, Cchar}([0 for i in 1:64]),
    0,0,0,
    ASI_FALSE,
    ASI_BAYER_RG,
    NTuple{16, Cint}([0 for i in 1:16]),
    NTuple{8, ASI_IMG_TYPE}([ASI_IMG_END for i in 1:8]),
    0,
    ASI_FALSE,
    ASI_FALSE,
    ASI_FALSE,
    ASI_FALSE,
    ASI_FALSE,
    1.,
    0,
    ASI_FALSE,
    NTuple{16, UInt8}([0 for i in 1:16])
)

# human-readable type
struct ASI_CAMERA_INFO
    Name::String
    CameraID::Cint
    MaxHeight::Clong
    MaxWidth::Clong
    IsColorCam::ASI_BOOL
    BayerPattern::ASI_BAYER_PATTERN
    SupportedBins::Vector
    SupportedVideoFormat::Vector
    PixelSize::Cdouble
    MechanicalShutter::ASI_BOOL
    ST4Port::ASI_BOOL
    IsCoolerCam::ASI_BOOL
    IsUSB3Host::ASI_BOOL
    IsUSB3Camera::ASI_BOOL
    ElecPerADU::Cfloat
    BitDepth::Cint
    IsTriggerCam::ASI_BOOL
end

ASI_CAMERA_INFO(info::_ASI_CAMERA_INFO) = ASI_CAMERA_INFO(
    split(String([Char(c) for c in info.Name]), "\0")[1],
    info.CameraID,
    info.MaxHeight,
    info.MaxWidth,
    info.IsColorCam,
    info.BayerPattern,
    [i for i in info.SupportedBins if i > 0],
    [i for i in info.SupportedVideoFormat if i > ASI_IMG_END],
    info.PixelSize,
    info.MechanicalShutter,
    info.ST4Port,
    info.IsCoolerCam,
    info.IsUSB3Host,
    info.IsUSB3Camera,
    info.ElecPerADU,
    info.BitDepth,
    info.IsTriggerCam
)

mutable struct _ASI_CONTROL_CAPS
    Name::NTuple{64, Cchar}
    Description::NTuple{128, UInt8}
    MaxValue::Clong
    MinValue::Clong
    DefaultValue::Clong
    IsAutoSupported::ASI_BOOL
    IsWritable::ASI_BOOL
    ControlType::ASI_CONTROL_TYPE
    Unused::NTuple{32, UInt8}
end

_ASI_CONTROL_CAPS() = _ASI_CONTROL_CAPS(
    NTuple{64, Cchar}([0 for i in 1:64]),   # Name
    NTuple{128, Cchar}([0 for i in 1:128]), # Description
    0, 0, 0,    # Max, Min, Default
    ASI_FALSE,  # IsAutoSupported
    ASI_FALSE,  # IsWritable
    ASI_GAIN,   # ControlType; gain corresponds to 0
    NTuple{32, UInt8}([0 for i in 1:32])    # Unused
)

# human-readable type
struct ASI_CONTROL_CAPS
    Name::String
    Description::String
    MaxValue::Clong
    MinValue::Clong
    DefaultValue::Clong
    IsAutoSupported::ASI_BOOL
    IsWritable::ASI_BOOL
    ControlType::ASI_CONTROL_TYPE
end

ASI_CONTROL_CAPS(control_caps::_ASI_CONTROL_CAPS) = ASI_CONTROL_CAPS(
    split(String([Char(c) for c in control_caps.Name]), "\0")[1],   # Name
    split(String([Char(c) for c in control_caps.Description]), "\0")[1], # Description
    control_caps.MaxValue, control_caps.MinValue, control_caps.DefaultValue,    # Max, Min, Default
    control_caps.IsAutoSupported,  # IsAutoSupported
    control_caps.IsWritable,       # IsWritable
    control_caps.ControlType       # ControlType; gain corresponds to 0
)

# const ASI_CONTROL_CAPS = _ASI_CONTROL_CAPS

mutable struct _ASI_ID
    id::NTuple{8, Cuchar}
end

_ASI_ID() = _ASI_ID(NTuple{8, Cuchar}(zeros(Cuchar, 10)))

const ASI_ID = _ASI_ID
const ASI_SN = ASI_ID

mutable struct _ASI_SUPPORTED_MODE
    SupportedCameraModes::NTuple{16, ASI_CAMERA_MODE}
end

_ASI_SUPPORTED_MODE() = _ASI_SUPPORTED_MODE(
    NTuple{16, ASI_CAMERA_MODE}(fill(ASI_MODE_END, 16))
)

const ASI_SUPPORTED_MODE = _ASI_SUPPORTED_MODE


text:
module JSort
export Parameters,
       linearalign, goodchannels, findpeaks, gainshift,
       leastsquares, localalign, splitregion, featurealign, boxbounds,
       arraybox, splitregions2d, featureoverlap2d, featurealign2d,
       alignspectra, bin, binhist, aligntime,
       histogram_edges, histogram, L1, L2, χ², Weight, ThresholdSum,
       ThresholdArea, featurealign2d_weighted, annotate_plot_2d, compare_2d,
    L1l, L2l, neighbourreduce, plotresiduals, compare_methods, splitmetric,
    ParticleDetector, ParticleCalibrator, read, readparticledetector,
    fill!, combine!, addcalibrator!, CBuffer, CBuffer2D, read, write!, next!, close,
    open, isgood, plot, calibrateall, invert, getindex

include("gainshift.jl")
export gainshift

include("./BetheBloch/BetheBloch.jl")
export Banana

include("./buffer/buffer.jl")
include("gate.jl")
export Gate, center, gate, center, Point, Gate2D, AbstractGate

include("cascade.jl")
export GammaCascade, extend!

include("levels.jl")
export fetchlevels

include("calibrator.jl")
export TimeFunction, TimeCalibrator, GammaCalibrator
       TimeCorrector, AbstractCalibrator, AbstractCalibrator1D

include("peakcombine.jl")
export combinepeaks

include("detectors.jl")
export GammaDetector, TimeDetector, readres, Detector, ChannelDetector,
       align, align!

include("parameters.jl")
include("metrics.jl")

include("weights.jl")
export ThresholdSumShift, WellSum

include("regression.jl")
include("reader.jl")
export AbstractReader, Reader, Event

include("quickreader.jl")
export QuickReader

include("alignment.jl")

include("alignment2d.jl")
export automap

include("regionalign.jl")
export regionpeaks, regionalign, backshift, plothist

include("iterativealign.jl")
export iterativealign, GammaCascade

include("histogram.jl")
include("timealignment.jl")

include("interactive.jl")
export GatePicker, pickpgate

#include("BetheBloch/BetheBloch.jl")
#using .BetheBloch
#export Material, Particle, Banana

end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using DataFrames: DataFrame
using Dates
using TimeZones

@testset "provide" begin
    @testset "basic" begin
        @system SProvide(Controller) begin
            a => DataFrame(index=(0:3)u"hr", value=0:10:30) ~ provide
        end
        s = instance(SProvide)
        @test s.a'.index == [0, 1, 2, 3]u"hr"
        @test s.a'.value == [0, 10, 20, 30]
    end

    @testset "index" begin
        @system SProvideIndex(Controller) begin
            a => DataFrame(i=(0:3)u"hr", value=0:10:30) ~ provide(index=:i)
        end
        s = instance(SProvideIndex)
        @test s.a'.i == [0, 1, 2, 3]u"hr"
    end

    @testset "autounit" begin
        @system SProvideAutoUnit(Controller) begin
            a => DataFrame("index (hr)" => 0:3, "value (m)" => 0:10:30) ~ provide
            b => DataFrame("index" => (0:3)u"hr", "value (m)" => 0:10:30) ~ provide(autounit=false)
        end
        s = instance(SProvideAutoUnit)
        @test s.a'."index" == [0, 1, 2, 3]u"hr"
        @test s.a'."value" == [0, 10, 20, 30]u"m"
        @test s.b'."index" == [0, 1, 2, 3]u"hr"
        @test s.b'."value (m)" == [0, 10, 20, 30]
    end

    @testset "time" begin
        @system SProvideTime(Controller) begin
            a => DataFrame("index (d)" => 0:3, "value (m)" => 0:10:30) ~ provide
        end
        c = :Clock => :step => 1u"d"
        s = instance(SProvideTime; config=c)
        @test s.a'.index == [0, 1, 2, 3]u"d"
        @test s.a'.value == [0, 10, 20, 30]u"m"
        @test_throws ErrorException instance(SProvideTime)
    end

    @testset "tick" begin
        @system SProvideTick(Controller) begin
            a => DataFrame(index=0:3, value=0:10:30) ~ provide(init=context.clock.tick, step=1)
        end
        s = instance(SProvideTick)
        @test s.a'.index == [0, 1, 2, 3]
    end

    @testset "calendar" begin
        @system SProvideCalendar(Controller) begin
            calendar(context) ~ ::Calendar
            a ~ provide(init=calendar.time, parameter)
        end
        t0 = ZonedDateTime(2011, 10, 29, 0, tz"Asia/Seoul")
        t1 = ZonedDateTime(2011, 10, 29, 1, tz"Asia/Seoul")
        t2 = ZonedDateTime(2011, 10, 29, 3, tz"Asia/Seoul")
        Δt = Hour(1)
        df = DataFrame(index=t0:Δt:t2, value=0:3)
        c = (:Calendar => :init => t1, :0 => :a => df)
        s = instance(SProvideCalendar; config=c)
        @test s.a'.index == t1:Δt:t2
        @test s.a'.value == 1:3
    end

    @testset "parameter" begin
        @system SProvideParameter(Controller) begin
            a ~ provide(parameter)
        end
        df = DataFrame("index (hr)" => 0:3, "value (m)" => 0:10:30)
        c = :0 => :a => df
        s = instance(SProvideParameter; config=c)
        @test s.a'.index == [0, 1, 2, 3]u"hr"
        @test s.a'.value == [0, 10, 20, 30]u"m"
    end

    @testset "csv" begin
        @system SProvideCSV(Controller) begin
            a ~ provide(parameter)
        end
        s = mktemp() do f, io
            write(io, "index (hr),value\n0,0\n1,10\n2,20\n3,30")
            close(io)
            c = :0 => :a => f
            instance(SProvideCSV; config=c)
        end
        @test s.a'.index == [0, 1, 2, 3]u"hr"
        @test s.a'.value == [0, 10, 20, 30]
    end
end


text:
function toh05()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
    k = [(-19:0)...]
    z = exp2.(k)
    p = poly(z)
    z = [z ones(20)]

    p, PolyZeros(z)
end


text:
using NearestNeighbors, StaticArrays, Statistics
using DelayEmbeddings
using DelayEmbeddings: AbstractDataset

export AbstractLocalModel
export AverageLocalModel,LinearLocalModel
export localmodel_tsp, localmodel_cp
export MSEp

"""
    AbstractLocalModel
Supertype of methods for making a prediction of a query point `q` using local models,
following the methods of [1]. Concrete subtypes are `AverageLocalModel` and
`LinearLocalModel`.

All models weight neighbors with a chosen function, so that distant neighbors
have smaller impact on the prediction and so that the interpolation
is smooth. The default weighting function we use is
```math
\\begin{aligned}
ω_i(d_i,d_{max}) = \\left[ 1- \\left(\\frac{d_i}{d_{max}}\\right)^2\\right]^4
\\end{aligned}
```
with ``d_i = ||x_{nn,i} -q||_2`` being the distance of each neighbor from
the query point.

You can also provide your own function or give
`ω_safe(d, dmax) = dmax > 0 ? (1.1 - (d/dmax)^2)^4 : 1.0`
for a safe version of ``ω`` that takes into acount edge cases. Finally
you can also give `nothing` in place of `ω`. In that case no weighting is done
and direct average of neighbors is returned.

### Average Local Model

    AverageLocalModel(ω)

The prediction is simply the weighted average of the images ``y_{nn, i}`` of
the neighbors ``x_{nn, i}`` of the query point `q`, weighting using given
function `ω`
```math
\\begin{aligned}
y_{pred} = \\frac{\\sum{\\omega_i y_{nn,i}}}{\\sum{\\omega_i}}
\\end{aligned}
```

### Linear Local Model

    LinearLocalModel([ω ], μ::Real=2.])
    LinearLocalModel([ω ], s_min::Real, s_max::Real)

The prediction is a weighted linear regression over the neighbors ``x_{nn, i}`` of
the query and their images ``y_{nn,i}`` as shown in [1].

Giving either `μ` or `s_min` and `s_max` determines which type of regularization
is applied.
  * `μ` : Ridge Regression
    ```math
    \\begin{aligned}
    f(\\sigma) = \\frac{\\sigma^2}{\\mu^2 + \\sigma^2}
    \\end{aligned}
    ```
  *  `s_min`, `s_max` : Soft Threshold
    ```math
    \\begin{aligned}
    f(\\sigma) = \\begin{cases} 0, & \\sigma < s_{min}\\\\
    \\left(1 - \\left( \\frac{s_{max}-\\sigma}{s_{max}-s_{min}}\\right)^2 \\right)^2,
    &s_{min} \\leq \\sigma \\leq s_{max} \\\\
    1, & \\sigma > s_{max}\\end{cases}
    \\end{aligned}
    ```

## References
[1] : D. Engster & U. Parlitz, *Handbook of Time Series Analysis* Ch. 1,
VCH-Wiley (2006)
"""
abstract type AbstractLocalModel end

ω_safe(d, dmax) = dmax > 0 ? (x=d/dmax; y= 1.1-x*x; y*y*y*y) : 1.
ω_unsafe(d, dmax) = (x=d/dmax; y= 1-x*x; y*y*y*y)

"""
    AverageLocalModel(ω::Function = ω_unsafe)
See [`AbstractLocalModel`](@ref).
"""
struct AverageLocalModel{F} <: AbstractLocalModel
    ω::F
end
AverageLocalModel() = AverageLocalModel(ω_unsafe)

function (M::AverageLocalModel)(q,xnn,ynn,dists)
    if length(ynn) > 1 && (dmax = maximum(dists)) > 0
        y_pred = zeros(typeof(ynn[1]))
        Ω = zero(typeof(dmax))
        for (y,d) in zip(ynn,dists)
            ω2 = M.ω.(d, dmax)
            Ω += ω2
            y_pred = y_pred .+ ω2*y
        end
        y_pred /= Ω
        return y_pred
    end
    return ynn[1]
end

(M::AverageLocalModel{Nothing})(q,xnn,ynn,dists) = sum(ynn)/length(ynn)

"""
    LinearLocalModel([ω ], μ::Real=2.])
    LinearLocalModel([ω ], s_min::Real, s_max::Real)
See [`AbstractLocalModel`](@ref).
"""
struct LinearLocalModel{Ω, F} <: AbstractLocalModel
    ω::Ω # weighting
    f::F # regularization
end

LinearLocalModel(μ::Real=2.) = LinearLocalModel(ω_unsafe, μ)

#Regularization functions
ridge_reg(σ, μ) = σ^2/(μ^2 + σ^2)
function mcnames_reg(σ,smin,smax)
    if σ < smin return 0
    elseif σ > smax return 1
    else return (1-((smax-σ)/(smax-smin))^2)^2
    end
end

LinearLocalModel(ω,μ::Real=2.) =
LinearLocalModel(ω,(σ) -> ridge_reg(σ, μ))

LinearLocalModel(ω, s_min::Real, s_max::Real) =
LinearLocalModel(ω, (σ) -> mcnames_reg(σ, s_min, s_max))

function (M::LinearLocalModel)(
    q,
    xnn::Vector{SVector{L,T}},
    ynn::Vector{TT},
    dists) where {L,T,TT}

    @assert length(ynn)>0 "No Nearest Neighbors given"
    k= length(xnn)
    #Weight Function
    dmax = maximum(dists)
    #Create Weight Matrix
    W = Diagonal([M.ω.(di,dmax) for di in dists])
    x_mean = mean(xnn)
    y_mean = mean(ynn)
    #Create X
    X = zeros(k,L)
    for i=1:k
        X[i,1:end] = xnn[i] .- x_mean
    end
    #Pseudo Inverse
    U,S,V = svd(W*X)

    #Regularization
    #D+1 Singular Values
    Sp = Diagonal([σ>0 ? M.f(σ)/σ : 0 for σ in S])
    Xw_inv = V*Sp*U'
    #The following code is meant for 1D ynn values
    #Repeat for all components
    y_pred = map(eachindex(ynn[1])) do i
        y = map(ynn->ynn[i], ynn)
        #Coefficient Vector
        ν = Xw_inv * W* y
        y_mean[i] + (q-x_mean)'* ν
    end

    return y_pred
end

#= # THE FOLLOWING IS THE POLYNOMIAL MODEL, WHICH WE DID NOT IMPLEMENT
# BECAUSE IT WAS NOT NEARLY AS GOOD AS THE OTHER TWO.
# SOMEONE MIGHT WANT TO TAKE A LOOK AT ANY POINT IN THE FUTURE
struct LocalPolynomialModel <: AbstractLocalModel
    n::Int #n=0,1,2,3 Exponent in weighting function
    m::Int #m=1,2,3,4 degree of Polynome
end
function (M::LocalPolynomialModel)(
    q,
    xnn::Vector{SVector{D,T}},
    ynn,
    dists) where {D,T}

    @assert length(ynn)>0 "No Nearest Neighbors given"
    y_pred = zeros(size(ynn[1]))
    k= length(xnn) #Can this be inferred? Nope, and probably not worth it.
    #Weight Function
    ω(r) = (1-r^M.n)^M.n
    dmax = maximum(dists)
    #Create Weight Matrix
    W = diagm([ω(di/dmax) for di in dists])
    x_mean = mean(xnn)
    y_mean = mean(ynn)

    #Create X
    X = zeros(k,M.m*D)
    for i=1:k, j=1:M.m
        X[i,1+(j-1)*D:j*D] = (xnn[i] - x_mean).^j
    end
    #Pseudo Inverse
    U,S,V = svd(W*X)
    #Regularization
    #μ = 0.8
    #f(σ) = σ^2/(μ^2+σ^2)

    #This could be a field "Regularization function"
    #in LocalPolynomialModel
    smin = 0.08
    smax = 1
    function f(σ)
        if σ < smin return 0
        elseif σ > smin return 1
        else return (1-((smax-σ)/(smax-smin))^2)^2
        end
    end
    Sp = diagm([σ>0 ? f(σ)/σ : 0 for σ in S])
    Xw_inv = V*Sp*U'

    #The following code is meant for 1D ynn values
    #Repeat for all components
    for i=1:length(ynn[1])
        y = map(ynn->ynn[i], ynn)
        #Coefficient Vector
        ν = Xw_inv * W* y
        #Query Matrix
        q_mat = T[]; sizehint!(q_mat, M.m*D)
        for j=1:M.m
            append!(q_mat,(q-x_mean).^j)
        end
        y_pred[i] = y_mean[i] + dot(q_mat,ν)
    end

    return y_pred
end
=#

function neighborhood_and_distances(point::AbstractVector,R::AbstractDataset, tree,
                      ntype::FixedMassNeighborhood, n::Int, w::Int = 1)
    idxs,dists = knn(tree, point, ntype.K, false, i -> abs(i-n) < w)
    return idxs,dists
end
function neighborhood_and_distances(point::AbstractVector,R, tree,
    ntype::FixedMassNeighborhood)
    idxs,dists = knn(tree, point, ntype.K, false)
    return idxs,dists
end

function neighborhood_and_distances(point::AbstractVector,R::AbstractDataset, tree,
                      ntype::FixedSizeNeighborhood, n::Int, w::Int = 1)
    idxs = inrange(tree, point, ntype.ε)
    filter!((el) -> abs(el - n) ≥ w, idxs)
    dists = [norm(x-point) for x in R[idxs]]   #Note: this is not an SVector!
    return idxs,dists
end
function neighborhood_and_distances(point::AbstractVector,R::AbstractDataset,
                      tree, ntype::FixedSizeNeighborhood)
    idxs = inrange(tree, point, ntype.ε)
    dists = [norm(x-point) for x in R[idxs]]   #Note: this is not an SVector!
    return idxs,dists
end

function _localmodel_tsp(R::AbstractDataset{D,T},
                        tree::KDTree,
                        q::SVector{D,T},
                        p::Int;
                        method::AbstractLocalModel = AverageLocalModel(),
                        ntype::AbstractNeighborhood  = FixedMassNeighborhood(2),
                        stepsize::Int = 1) where {D,T}

    s_pred = Vector{SVector{D, T}}(undef, p+1)
    s_pred[1] = q
    for n=2:p+1   #Iteratively estimate timeseries
        idxs,dists = neighborhood_and_distances(q,R, tree,ntype)
        xnn = R[idxs]
        ynn = R[idxs.+stepsize]
        q = method(q, xnn, ynn, dists)
        s_pred[n] = q
    end
    return Dataset(s_pred)
end

"""
    localmodel_tsp(s, γ::Int, τ, p::Int; method, ntype, stepsize)
    localmodel_tsp(s, p::Int; method, ntype, stepsize)

Perform a timeseries prediction for `p` points,
using local weighted modeling [1]. The function always returns an
object of the same type as `s`, which can be either a timeseries (vector) or an
`AbstractDataset` (trajectory), and the returned data
always contains the final point of `s` as starting point. This means that the
returned data has length of `p + 1`.

If given `(s, γ, τ)`, it first calls `reconstruct` (from `DelayEmbeddings`) on `s`
with `(γ, τ)`. If given only `s` then no reconstruction is done.

## Keyword Arguments
  * `method = AverageLocalModel(ω_unsafe)` : Subtype of [`AbstractLocalModel`](@ref).
  * `ntype = FixedMassNeighborhood(2)` : Subtype of `AbstractNeighborhood` (from
    `DelayEmbeddings`).
  * `stepsize = 1` : Prediction step size.

## Description
Given a query point, the function finds its neighbors using neighborhood `ntype`.
Then, the neighbors `xnn` and their images `ynn` are used to make a prediction for
the future of the query point, using the provided `method`.
The images `ynn` are the points `xnn` shifted by `stepsize` into the future.

The algorithm is applied iteratively until a prediction of length `p` has
been created, starting with the query point to be the last point of the timeseries.

## References
[1] : D. Engster & U. Parlitz, *Handbook of Time Series Analysis* Ch. 1,
VCH-Wiley (2006)
"""
function localmodel_tsp(R::AbstractDataset{B}, p::Int;
    method::AbstractLocalModel = AverageLocalModel(),
    ntype::AbstractNeighborhood = FixedMassNeighborhood(2),
    stepsize::Int = 1) where B
    B > 1 || throw(ArgumentError("Dataset Dimension needs to be >1! ",
    "Alternatively pass embedding parameters."))
    return _localmodel_tsp(R, KDTree(R[1:end-stepsize]), R[end], p;
    method=method, ntype=ntype, stepsize=stepsize)
end

function localmodel_tsp(
    s::AbstractVector, γ::Int, τ::T, p::Int; kwargs... ) where {T}
    localmodel_tsp(reconstruct(s, γ, τ), p; kwargs...)[:,γ+1]
end

function localmodel_tsp(
    ss::AbstractDataset{B}, γ::Int, τ::T, p::Int; kwargs...) where {B,T}
    sind = SVector{B, Int}(((γ+1)*B - i for i in B-1:-1:0)...)
    localmodel_tsp(reconstruct(ss, γ, τ), p; kwargs...)[:,sind]
end

#####################################################################################
#                                  Cross Prediction                                   #
#####################################################################################

"""
    localmodel_cp(source_pool, target_pool, source_pred,  γ, τ; kwargs...)

Perform a cross prediction from  _source_ to _target_,
using local weighted modeling [1]. `source_pred` is the input for the prediction
and `source_pool` and `target_pool` are used as pooling/training data for the predictions.
The function always returns an object of the same type as `target_pool`,
which can be either a timeseries (vector) or an `AbstractDataset` (trajectory).

## Keyword Arguments
  * `method = AverageLocalModel(ω_unsafe)` : Subtype of [`AbstractLocalModel`](@ref).
  * `ntype = FixedMassNeighborhood(2)` : Subtype of `AbstractNeighborhood` (from
    `DelayEmbeddings`).
  * `stepsize = 1` : Prediction step size.

Instead of passing `γ` & `τ` for reconstruction one may also give
existing `Dataset`s as `source_pool` and `source_pred`.
In this case an additional keyword argument `y_idx_shift::Int=0` may be necessary
to account for the index shift introduced in the reconstruction process.

## Description
Given a query point, the function finds its neighbors using neighborhood `ntype`.
Then, the neighbors `xnn` and their images `ynn` are used to make a prediction for
the image of the query point, using the provided `method`.

## References
[1] : D. Engster & U. Parlitz, *Handbook of Time Series Analysis* Ch. 1,
VCH-Wiley (2006)
"""
function localmodel_cp(R::AbstractDataset{D,T},
                       target_train,
                       source_pred::AbstractDataset{D,T},
                       tree::KDTree;
                       method::AbstractLocalModel = AverageLocalModel(),
                       ntype::AbstractNeighborhood  = FixedMassNeighborhood(2),
                       y_idx_shift::Int=0) where {D,T}

    N = length(source_pred)
    target_pred = typeof(target_train)(undef, N)
    for n=1:N   #Iteratively estimate timeseries
        q = source_pred[n]
        idxs,dists = neighborhood_and_distances(q,R, tree,ntype)
        xnn = R[idxs]
        ynn = target_train[idxs .+ y_idx_shift]
        target_pred[n] = method(q, xnn, ynn, dists)[1]
    end
    return target_pred
end

function localmodel_cp(
    source_train::AbstractDataset{B},
    target_train,
    source_pred::AbstractDataset{B};kwargs...) where B
    B > 1 || throw(ArgumentError("Dataset Dimension needs to be >1! ",
    "Alternatively pass embedding parameters."))
    return localmodel_cp(source_train, target_train, source_pred,KDTree(source_train); kwargs... )
end

function localmodel_cp(
    source_train,
    target_train,
    source_pred,
     γ::Int, τ::Int; kwargs... )
    localmodel_cp(reconstruct(source_train, γ, τ),
                    target_train,
                    reconstruct(source_pred, γ, τ);
                    y_idx_shift=γ*τ, kwargs...)
end
function localmodel_cp(
    source_train,
    target_train,
    source_pred,
     γ::Int, τ::T; kwargs... ) where {T}
    localmodel_cp(reconstruct(source_train, γ, τ),
                    target_train,
                    reconstruct(source_pred, γ, τ);
                    y_idx_shift=maximum(τ), kwargs...)
end

#####################################################################################
#                                  Error Measures                                   #
#####################################################################################

function MSE1(
    R::AbstractDataset{D,T},
    tree::KDTree,
    R_test::AbstractDataset{D,T};
    method::AbstractLocalModel = AverageLocalModel(),
    ntype::AbstractNeighborhood  = FixedMassNeighborhood(2),
    stepsize::Int = 1) where {D,T}

    y_test = map(q-> q[end], R_test[2:end])
    y_pred = T[]; sizehint!(y_pred, length(y_test))
    for q in R_test[1:end-1] # Remove last state,
                            #because there is nothing to compare the preciction to
        idxs,dists = neighborhood_and_distances(q, R, tree,ntype)
        xnn = R[idxs]
        ynn = R[idxs.+stepsize]
        push!(y_pred,method(q,xnn,ynn, dists)[end])
    end
    return norm(y_test-y_pred)^2/length(y_test)
end
#FIXME: I shouldn't have to square the norm... What is the solution?

"""
    MSEp(R::AbstractDataset{D,T}, R_test, p; method, ntype, stepsize) -> error

Compute mean squared error of iterated predictions of length `p`
using test set `R_test`.

## Description
This error measure takes in a prediction model consisting of `R`,
`method`, `ntype` and `stepsize` and evaluates its performance. The test set
`R_test` is
a delay reconstruction with the same delay `τ` and dimension `D` as `R`.
For each subset of `R_test` with length `p` it calls `localmodel_tsp`.
The model error is then defined as
```math
\\begin{aligned}
MSE_p = \\frac{1}{p|T_{ref}|}\\sum_{t\\in T_{ref}}\\sum_{i=1}^{p} \\left(y_{t+i}
- y_{pred,t+i} \\right)^2
\\end{aligned}
```
where ``|T_{ref}|`` is the number of subsets of `R_test` used.
## References
See [`localmodel_tsp`](@ref).
"""
function MSEp(
    R::AbstractDataset{D,T},
    tree::KDTree,
    R_test::AbstractDataset{D,T},
    p::Int;
    kwargs...) where {D,T}
    if p == 1
        return MSE1(R,tree,R_test;kwargs...)
    end
    Tref = (length(R_test)-p-1)
    error = 0
    for t =1:Tref
        R_pred = _localmodel_tsp(R,tree,R_test[t], p; kwargs...)
        error += norm(R_test[t:t+p]-R_pred.data)^2 /Tref/p
    end
    return error
end
#FIXME: I shouldn't have to square the norm... What is the solution?
MSEp(R, R_test, p; method::AbstractLocalModel = AverageLocalModel(),
ntype::AbstractNeighborhood = FixedMassNeighborhood(2),
stepsize::Int = 1) =
MSEp(R, KDTree(R[1:end-stepsize]), R_test, p; method=method,
ntype=ntype, stepsize=stepsize)


text:
@testset "scanmle KS statistic" begin
    seed = 11
    α = 0.5
    data = Example.makeparetodata(α, seed)
    @test MLE.scanKS(data, range(.4, length=11, stop=.6)) == [0.48, 0.5, 0.52]
    @test abs(MLE.scanmle(data).alpha - 1.5) < 1e-2
end


text:
# Author: Ashwani Rathee
# Program to Calculate Sum of Natural Numbers

println("Enter n: ")

n=readline()
n=parse(Float64,n)
sum=(n^2+n)/2
println("The sum of n-natural terms is ",sum)

max_stars_repo_name: ashwani-rathee/Learning--Julia

max_issues_repo_name: ashwani-rathee/Learning--Julia

max_forks_repo_name: ashwani-rathee/Learning--Julia


text:
using FastTransforms, Compat
using Compat.Test, Compat.LinearAlgebra

@testset "Chebyshev--Legendre transform" begin
    for k in round.([Int],10 .^ range(1,stop=4,length=20))
        r = randn(k)./(√).(1:k) # Proven 𝒪(√(log N)) error for ASY method.
        @test leg2cheb(r) ≈ cjt(r,0.,0.)
    end

    @test leg2chebu([1.0,2,3,4,5])  ≈ [0.546875,0.5,0.5390625,1.25,1.3671875]

    c = randn(1000)./(√).(1:1000);

    @test leg2cheb(cheb2leg(c)) ≈ c
    @test cheb2leg(leg2cheb(c)) ≈ c

    @test norm(jac2jac(c,0.,√2/2,-1/4,√2/2)-jjt(c,0.,√2/2,-1/4,√2/2),Inf) < 10length(c)*eps()

    @test norm(ultra2ultra(ultra2ultra(c,.5,.75),.75,.5)-c,Inf) < 10length(c)*eps()

    @test FastTransforms.th_cheb2leg([1.0,2,3,4,5]) ≈ cheb2leg([1.0,2,3,4,5])
    @test FastTransforms.th_leg2cheb([1.0,2,3,4,5]) ≈ leg2cheb([1.0,2,3,4,5])
end


text:
using Test
using StaticArrays
using VariableTemplates

@testset "Complex models" begin

    #### Test models
    #
    #   EmptyModel
    #   ScalarModel
    #   VectorModel
    #   MatrixModel
    #   CompositModel
    #   NTupleContainingModel

    abstract type OneLayerModel end

    struct EmptyModel <: OneLayerModel end
    vars_state(m::EmptyModel, T) = @vars()

    struct ScalarModel <: OneLayerModel end
    vars_state(m::ScalarModel, T) = @vars(x::T)

    struct VectorModel{N} <: OneLayerModel end
    vars_state(m::VectorModel{N}, T) where {N} = @vars(x::SVector{N, T})

    struct MatrixModel{N,M} <: OneLayerModel end
    vars_state(m::MatrixModel{N,M}, T) where {N,M} = @vars(x::SHermitianCompact{N, T, M})

    abstract type TwoLayerModel end

    Base.@kwdef struct CompositModel{Nv,N,M} <: TwoLayerModel
        empty_model = EmptyModel()
        scalar_model = ScalarModel()
        vector_model = VectorModel{Nv}()
        matrix_model = MatrixModel{N,M}()
    end
    function vars_state(m::CompositModel, T)
        @vars begin
            empty_model::vars_state(m.empty_model, T)
            scalar_model::vars_state(m.scalar_model, T)
            vector_model::vars_state(m.vector_model, T)
            matrix_model::vars_state(m.matrix_model, T)
        end
    end

    Base.@kwdef struct NTupleModel <: OneLayerModel
        scalar_model = ScalarModel()
    end

    function vars_state(m::NTupleModel, T)
        @vars begin
            scalar_model::vars_state(m.scalar_model, T)
        end
    end

    vars_state(m::NTuple{N, NTupleModel}, FT) where {N} =
        Tuple{ntuple(i -> vars_state(m[i], FT), N)...}

    Base.@kwdef struct NTupleContainingModel{N,Nv} <: TwoLayerModel
        ntuple_model = ntuple(i -> NTupleModel(), N)
        vector_model = VectorModel{Nv}()
        scalar_model = ScalarModel()
    end

    function vars_state(m::NTupleContainingModel, T)
        @vars begin
            ntuple_model::vars_state(m.ntuple_model, T)
            vector_model::vars_state(m.vector_model, T)
            scalar_model::vars_state(m.scalar_model, T)
        end
    end

    FT = Float32;

    # ------------------------------- Test concretize
    m = EmptyModel()
    st = vars_state(m, FT)
    nt = concretize(st)
    @test nt == NamedTuple()

    m = ScalarModel()
    st = vars_state(m, FT)
    nt = concretize(st)
    @test nt.x.type == Val(FT)
    @test nt.x.gid == Val(1)
    @test varsize(nt) == 1

    Nv = 3
    m = VectorModel{Nv}()
    st = vars_state(m, FT)
    nt = concretize(st)
    @test nt.x.type == Val(SVector{Nv,FT})
    @test nt.x.gid == Val(1)
    @test varsize(nt) == 3

    N = 3
    M = 4
    m = MatrixModel{N,M}()
    st = vars_state(m, FT)
    nt = concretize(st)
    @test nt.x.type == Val(SHermitianCompact{N, FT, M})
    @test nt.x.gid == Val(1)
    @test varsize(nt) == M

    Nv = 3
    N = 3
    M = 4
    m = CompositModel{Nv,N,M}()
    st = vars_state(m, FT)
    nt = concretize(st)
    @test nt.empty_model == NamedTuple()
    @test nt.scalar_model.x.gid == Val(1)
    @test nt.scalar_model.x.type == Val(FT)
    @test nt.vector_model.x.gid == Val(2)
    @test nt.vector_model.x.type == Val(SVector{Nv, FT})
    @test nt.matrix_model.x.gid == Val(5)
    @test nt.matrix_model.x.type == Val(SHermitianCompact{N, FT, M})
    @test varsize(nt) == M+1+Nv

    Nv = 3
    N = 3
    m = NTupleContainingModel{N,Nv}()
    st = vars_state(m, FT)
    for i in 1:N
        @test m.ntuple_model[i] isa NTupleModel
        @test m.ntuple_model[i].scalar_model isa ScalarModel
    end
    nt = concretize(st)
    @test length(nt.ntuple_model) == N
    for i in 1:N
        @test hasproperty(nt.ntuple_model[i], :scalar_model)
    end

    # ------------------------------- Test getproperty
    m = ScalarModel()
    st = vars_state(m, FT)
    nt = concretize(st)
    vs = varsize(nt)
    a_global = collect(1:vs)
    v = Vars{nt}(a_global)
    @test v.x == FT(1)

    Nv = 3
    m = VectorModel{Nv}()
    st = vars_state(m, FT)
    nt = concretize(st)
    vs = varsize(nt)
    a_global = collect(1:vs)
    v = Vars{nt}(a_global)
    @test v.x == SVector{Nv,FT}(FT[1,2,3])

    N = 3
    M = 6
    m = MatrixModel{N,M}()
    st = vars_state(m, FT)
    nt = concretize(st)
    vs = varsize(nt)
    a_global = collect(1:vs)
    v = Vars{nt}(a_global)
    @test v.x == SHermitianCompact{N,FT,M}(collect(1:1+M-1))

    Nv = 3
    N = 3
    M = 6
    m = CompositModel{Nv,N,M}()
    st = vars_state(m, FT)
    nt = concretize(st)
    vs = varsize(nt)
    a_global = collect(1:vs)
    v = Vars{nt}(a_global)

    scalar_model = v.scalar_model
    @test VariableTemplates.get_tup_chain(scalar_model) == (:scalar_model,)
    @test v.scalar_model.x == FT(1)

    vector_model = v.vector_model
    @test VariableTemplates.get_tup_chain(vector_model) == (:vector_model,)
    @test v.vector_model.x == SVector{Nv,FT}([2,3,4])

    matrix_model = v.matrix_model
    @test VariableTemplates.get_tup_chain(matrix_model) == (:matrix_model,)
    @test v.matrix_model.x == SHermitianCompact{N,FT,M}(collect(5:5+M-1))

    Nv = 3
    N = 3
    m = NTupleContainingModel{N,Nv}()
    st = vars_state(m, FT)
    nt = concretize(st)
    vs = varsize(nt)
    a_global = collect(1:vs)
    @test length(nt.ntuple_model) == N
    for i in 1:N
        @test hasproperty(nt.ntuple_model[i], :scalar_model)
        @test hasproperty(nt.ntuple_model[i], :scalar_model)
    end
    v = Vars{nt}(a_global)

    @test VariableTemplates.get_tup_chain(v.vector_model) == (:vector_model,)
    @test v.vector_model.x == SVector{Nv,FT}([4,5,6])

    @test VariableTemplates.get_tup_chain(v.scalar_model) == (:scalar_model,)
    @test v.scalar_model.x == FT(7)

    @test VariableTemplates.get_tup_chain(v.ntuple_model) == (:ntuple_model,)
    for i in 1:N
        @test VariableTemplates.get_tup_chain(v.ntuple_model[i]) == (:ntuple_model, i)
    end
    for i in 1:N
        @test VariableTemplates.get_tup_chain(v.ntuple_model[i].scalar_model) == (:ntuple_model, i, :scalar_model)
    end

    for i in 1:N
        @test v.ntuple_model[i].scalar_model.x == FT(i)
        @test v.vector_model.x == SVector{Nv,FT}(N+1:N+Nv)
        @test v.scalar_model.x == FT(N+Nv+1)
    end

    fn = flattenednames(st)
    @test fn[1] === "ntuple_model[1].scalar_model.x"
    @test fn[2] === "ntuple_model[2].scalar_model.x"
    @test fn[3] === "ntuple_model[3].scalar_model.x"
    @test fn[4] === "vector_model.x[1]"
    @test fn[5] === "vector_model.x[2]"
    @test fn[6] === "vector_model.x[3]"
    @test fn[7] === "scalar_model.x"

    ftc = flattened_tup_chain(st)
    @test ftc[1] === (:ntuple_model, 1, :scalar_model, :x)
    @test ftc[2] === (:ntuple_model, 2, :scalar_model, :x)
    @test ftc[3] === (:ntuple_model, 3, :scalar_model, :x)
    @test ftc[4] === (:vector_model, :x)
    @test ftc[5] === (:scalar_model, :x)

    # getproperty with tup-chain
    for i in 1:N
        @test v.ntuple_model[i].scalar_model.x == getproperty(v, (:ntuple_model, i, :scalar_model, :x))
        @test v.vector_model.x == getproperty(v, (:vector_model, :x))
        @test v.scalar_model.x == getproperty(v, (:scalar_model, :x))
    end

end


text:
"""
    readfasta(file::String);

Read a disk file of fasta sequence and return a named tuple of sequences.
Sequence names will be the first field of the id line, which is space, tab, comma, and pipe delimited. For instance: >2 sample2|Jun4,2021 will be named "2".
(Remember to end with ; when reading large files \U1F601)
"""
function readfasta(file::String)
    open(file) do f

        a = [i for i in eachline(f)]
        chr = Array{Symbol,1}()
        seq = Array{String,1}()

        for i in 1:length(a) - 1

            if i % 2 == 0
                continue
            end

            if !occursin(r"^>", a[i])
                error("Fasta file is not properly formatted: Line $i")
            end

            c = split(a[i], r">|\s+|\||,|\t")

            if length(c) < 2
                continue
            end

            c = Symbol(c[2])
            push!(chr, c)
            push!(seq, a[i+1])

        end

        chr = Tuple(chr)
        fs = NamedTuple{chr}(seq)

        return fs

    end
end


text:
# Load packages and helper functions
using MPI, Statistics
include("../parallelhelpers.jl")

# Call MPI.init() and
MPI.Init()
comm = MPI.COMM_WORLD
my_rank = MPI.Comm_rank(comm)
n_procs = MPI.Comm_size(comm)
root = 0

# Load the file names
dataLoc = "../../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Distribute filenames by rank
myfnames = fnames[my_rank+1:n_procs:length(fnames)]

# Calculate the word counts
mycounts =  getcounts.(myfnames)

# Gather all counts onto Process 0

# First we have to serialize the dictionary so we can send it
msg = MPI.serialize(mycounts)

# Since each msg size is different we have to get the length of each one
counts = Int32.(MPI.Allgather(length(msg),comm))

# Then we can gather on process 0 (root)
res = MPI.Gatherv(msg, counts, root, comm)

# And deserialize and merge the counts to get the overall word counts
if my_rank == 0
    idx = [0; cumsum(counts)]
    res = [res[idx[i]+1:idx[i+1]] for i in 1:length(idx)-1]
    allcounts = vcat(MPI.deserialize.(res)...)
    overallcounts = merge(+,allcounts...)
else
    overallcounts = nothing
end

# Then we can use Broadcast to send these counts to the rest of the processes
# (note the function bcast will do the serialization/deserialization for you, gather does not currently have an equivalent)
overallcounts = MPI.bcast(overallcounts, root, comm)

# Calculate the normalized counts
normcount = getnormcount.(mycounts,Ref(overallcounts))

# Print out the top 5 words for each file and their counts on the leader process
for i=1:length(normcount)
    top5 = sort(collect(normcount[i]), by=x->x[2], rev=true)[1:min(5,length(normcount[i])),:]
    println(top5)
end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
# --- Jacobian Types --- #

"""
    AbstractJacobian

Supertype for gradient/jacobian wrappers.
"""
abstract type AbstractJacobian end

"""
    Empty <: AbstractJacobian

A gradient/jacobian wrapper which indicates that a given gradient/jacobian is empty.
"""
struct Empty <: AbstractJacobian end

"""
    Zeros <: AbstractJacobian

A gradient/jacobian wrapper which indicates that a given gradient/jacobian is zero.
"""
struct Zeros <: AbstractJacobian end

"""
    Identity <: AbstractJacobian

A gradient/jacobian wrapper which indicates that a given gradient/jacobian is the identity
matrix.
"""
struct Identity <: AbstractJacobian end

"""
    Invariant <: AbstractJacobian

A gradient/jacobian wrapper which indicates that the value of given gradient/jacobian is
independent of the state rates, states, inputs, parameters, and time.

The value of the invariant jacobian/gradient may be provided upon construction, otherwise
it will be calculated using automatic differentiation.
"""
struct Invariant{T} <: AbstractJacobian; value::T; end

Invariant() = Invariant(nothing)

"""
    Constant <: AbstractJacobian

A gradient/jacobian wrapper which indicates that the value of given gradient/jacobian is
independent of the state rates, states, inputs, and time.

The jacobian/gradient will be calculated using automatic differentiation.  Alternatively,
a function may be provided which defines the value of the gradient/jacobian as a function
of the parameters.
"""
struct Constant{T} <: AbstractJacobian; func::T; end

Constant() = Constant(nothing)

"""
    Linear <: AbstractJacobian

A gradient/jacobian wrapper which indicates that the value of given gradient/jacobian is
linear with respect to the differentiated variable.

The jacobian/gradient will be calculated using automatic differentiation.  Alternatively,
a function may be provided which defines the value of the gradient/jacobian as a function
of the state rates, states, inputs, parameters, and time.
"""
struct Linear{T} <: AbstractJacobian; func::T; end

Linear() = Linear(nothing)

"""
    Nonlinear <: AbstractJacobian

A gradient/jacobian wrapper which indicates that the value of given gradient/jacobian is
nonlinear with respect to the differentiated variable.

The jacobian/gradient will be calculated using automatic differentiation.  Alternatively,
a function may be provided which defines the value of the gradient/jacobian as a function
of the state rates, states, inputs, parameters, and time.
"""
struct Nonlinear{T} <: AbstractJacobian; func::T; end

Nonlinear() = Nonlinear(nothing)

# trait dispatch
isempty(::AbstractJacobian) = false
isempty(::Empty) = true

iszero(::AbstractJacobian) = false
iszero(::Zeros) = true

isidentity(::AbstractJacobian) = false
isidentity(::Identity) = true

isinvariant(::AbstractJacobian) = false
isinvariant(::Empty) = true
isinvariant(::Zeros) = true
isinvariant(::Identity) = true
isinvariant(::Invariant) = true

isconstant(::AbstractJacobian) = false
isconstant(::Empty) = true
isconstant(::Zeros) = true
isconstant(::Identity) = true
isconstant(::Invariant) = true
isconstant(::Constant) = true

islinear(::AbstractJacobian) = false
islinear(::Empty) = true
islinear(::Zeros) = true
islinear(::Identity) = true
islinear(::Invariant) = true
islinear(::Constant) = true
islinear(::Linear) = true

# --- Jacobian Update Functions --- #

for iarg = (1, 2, 3, 4, 5)

    if iarg == 1
        fname = :add_rate_jacobian
        adfunc = :autodiff_jacobian_func
        ncol = :nx
    elseif iarg == 2
        fname = :add_state_jacobian
        adfunc = :autodiff_jacobian_func
        ncol = :nx
    elseif iarg == 3
        fname = :add_input_jacobian
        adfunc = :autodiff_jacobian_func
        ncol = :ny
    elseif iarg == 4
        fname = :add_parameter_jacobian
        adfunc = :autodiff_jacobian_func
        ncol = :np
    elseif iarg == 5
        fname = :add_time_gradient
        adfunc = :autodiff_derivative_func
        ncol = 1
    end

    @eval begin
        # by default, return original rate jacobian
        $(fname)(jac, f, iip, nx, ny, np) = jac

        # invariant jacobian
        function $(fname)(::Invariant{Nothing}, f, iip, nx, ny, np)

            # set values arbitrarily
            dx = zeros(val(nx))
            x = zeros(val(nx))
            y = zeros(val(ny))
            p = zeros(val(np))
            t = 0.0

            if iip
                # initialize jacobian
                J = zeros(nx, $ncol)
                # get in-place function
                fJ = $(adfunc)(f, val(nx), $iarg)
                # calculate in-place jacobian
                fJ(J, dx, x, y, p, t)
            else # out-of-place
                # get out-of-place function
                fJ = $(adfunc)(f, $iarg)
                # calculate out-of-place jacobian
                J = fJ(dx, x, y, p, t)
            end

            return Invariant(J)
        end

        # constant jacobian
        function $(fname)(::Constant{Nothing}, f, iip, nx, ny, np)

            # set values arbitrarily
            dx = zeros(nx)
            x = zeros(nx)
            y = zeros(ny)
            t = 0.0

            if iip
                # get in-place function
                fJ = $(adfunc)(f, val(nx), $iarg)
                # get in-place jacobian function
                fJ = (J, p) -> fJ(J, dx, x, y, p, t)
            else
                # get out-of-place function
                fJ = $(adfunc)(f, $iarg)
                # get out-of-place jacobian function
                fJ = (p) -> fJ(dx, x, y, p, t)

            end

            return Constant(fJ)
        end

        # out-of-place or in-place linear rate jacobian
        function $(fname)(::Linear{Nothing}, f, iip, nx, ny, np)
            if iip
                fJ = $(adfunc)(f, nx, $iarg)
            else
                fJ = $(adfunc)(f, $iarg)
            end
            return Linear(fJ)
        end

        # out-of-place or in-place nonlinear rate jacobian
        function $(fname)(::Nonlinear{Nothing}, f, iip, nx, ny, np)
            if iip
                fJ = $(adfunc)(f, nx, $iarg)
            else
                fJ = $(adfunc)(f, $iarg)
            end
            return Nonlinear(fJ)
        end

    end
end

for iarg = (1, 2, 3, 4)

    if iarg == 1
        fname = :add_coupling_rate_jacobian
        adfunc = :autodiff_jacobian_func
        ncol = :nx
    elseif iarg == 2
        fname = :add_coupling_state_jacobian
        adfunc = :autodiff_jacobian_func
        ncol = :nx
    elseif iarg == 3
        fname = :add_coupling_parameter_jacobian
        adfunc = :autodiff_jacobian_func
        ncol = :np
    elseif iarg == 4
        fname = :add_coupling_time_gradient
        adfunc = :autodiff_derivative_func
        ncol = 1
    end

    @eval begin
        # by default, return original rate jacobian
        $(fname)(jac, f, iip, nx, ny, np) = jac

        # invariant jacobian
        function $(fname)(::Invariant{Nothing}, f, iip, nx, ny, np)

            # set values arbitrarily
            dx = zeros(val(nx))
            x = zeros(val(nx))
            p = zeros(val(np))
            t = 0.0

            if iip
                # initialize jacobian
                J = zeros(nx, $ncol)
                # get in-place function
                fJ = $(adfunc)(f, val(ny), $iarg)
                # calculate in-place jacobian
                fJ(J, dx, x, p, t)
            else # out-of-place
                # get out-of-place function
                fJ = $(adfunc)(f, $iarg)
                # calculate out-of-place jacobian
                J = fJ(dx, x, p, t)
            end

            return Invariant(J)
        end

        # constant jacobian
        function $(fname)(::Constant{Nothing}, f, iip, nx, ny, np)

            # set values arbitrarily
            dx = zeros(nx)
            x = zeros(nx)
            t = 0.0

            if iip
                # get in-place function
                fJ = $(adfunc)(f, val(ny), $iarg)
                # get in-place jacobian function
                fJ = (J, p) -> fJ(J, dx, x, p, t)
            else
                # get out-of-place function
                fJ = $(adfunc)(f, $iarg)
                # get out-of-place jacobian function
                fJ = (p) -> fJ(dx, x, p, t)

            end

            return Constant(fJ)
        end

        # out-of-place or in-place linear rate jacobian
        function $(fname)(::Linear{Nothing}, f, iip, nx, ny, np)
            if iip
                fJ = $(adfunc)(f, val(ny), $iarg)
            else
                fJ = $(adfunc)(f, $iarg)
            end
            return Linear(fJ)
        end

        # out-of-place or in-place nonlinear rate jacobian
        function $(fname)(::Nonlinear{Nothing}, f, iip, nx, ny, np)
            if iip
                fJ = $(adfunc)(f, val(ny), $iarg)
            else
                fJ = $(adfunc)(f, $iarg)
            end
            return Nonlinear(fJ)
        end

    end
end

# --- Jacobian Evaluation --- #

# out-of-place, statically sized

function jacobian(::Empty, ::Val{NX}, ::Val{NY}, args...) where {NX, NY}
    return SMatrix{NX, NY, Float64}()
end

function jacobian(::Zeros, ::Val{NX}, ::Val{NY}, args...) where {NX, NY}
    return zeros(SMatrix{NX, NY, Float64})
end

function jacobian(::Identity, ::Val{NX}, ::Val{NY}, args...) where {NX, NY}
    return SMatrix{NX, NY}(I)
end

function jacobian(jac::Invariant, ::Val{NX}, ::Val{NY}, args...) where {NX, NY}
    return SMatrix{NX, NY}(jac.value)
end

function jacobian(jac::Constant, ::Val{NX}, ::Val{NY}, p) where {NX, NY}
    return SMatrix{NX, NY}(jac.func(p))
end

function jacobian(jac::Constant, ::Val{NX}, ::Val{NY}, dx, x, p, t) where {NX, NY}
    return SMatrix{NX, NY}(jac.func(p))
end

function jacobian(jac::Constant, ::Val{NX}, ::Val{NY}, dx, x, y, p, t) where {NX, NY}
    return SMatrix{NX, NY}(jac.func(p))
end

function jacobian(jac::Union{Linear, Nonlinear}, ::Val{NX}, ::Val{NY}, args...) where {NX, NY}
    return SMatrix{NX, NY}(jac.func(args...))
end

# out-of-place, dynamically sized

function jacobian(::Empty, nx, ny, args...)
    return Matrix{Float64}(nx, ny)
end

function jacobian(::Zeros, nx, ny, args...)
    return zeros(nx, ny)
end

function jacobian(::Identity, nx, ny, args...)
    return I(min(nx, ny))
end

function jacobian(jac::Invariant, nx, ny, args...)
    return jac.value
end

function jacobian(jac::Constant, nx, ny, p)
    return jac.func(p)
end

function jacobian(jac::Constant, nx, ny, dx, x, p, t)
    return jac.func(p)
end

function jacobian(jac::Constant, nx, ny, dx, x, y, p, t)
    return jac.func(p)
end

function jacobian(jac::Union{Linear, Nonlinear}, nx, ny, args...)
    return jac.func(args...)
end

# in-place

function jacobian!(J, ::Empty, args...)
    return J
end

function jacobian!(J, ::Zeros, args...)
    return J .= 0
end

function jacobian!(J, ::Identity, args...)
    J .= 0
    for i = 1:size(J,1)
        J[i,i] = 1
    end
    return J
end

function jacobian!(J, jac::Invariant, args...)
    return J .= jac.value
end

function jacobian!(J, jac::Constant, p)
    return jac.func(J, p)
end

function jacobian!(J, jac::Constant, dx, x, p, t)
    return jac.func(J, p)
end

function jacobian!(J, jac::Constant, dx, x, y, p, t)
    return jac.func(J, p)
end

function jacobian!(J, jac::Union{Linear, Nonlinear}, args...)
    return jac.func(J, args...)
end

# --- Gradient Evaluation --- #

# out-of-place, statically sized

function gradient(::Empty, ::Val{NX}, args...) where NX
    return SVector{0, Float64}()
end

function gradient(::Zeros, ::Val{NX}, args...) where NX
    return zeros(SVector{NX, Float64})
end

function gradient(jac::Invariant, ::Val{NX}, args...) where NX
    return SVector{NX}(jac.value)
end

function gradient(jac::Constant, ::Val{NX}, p) where NX
    return SVector{NX}(jac.func(p))
end

function gradient(jac::Constant, ::Val{NX}, dx, x, p, t) where NX
    return SVector{NX}(jac.func(p))
end

function gradient(jac::Constant, ::Val{NX}, dx, x, y, p, t) where NX
    return SVector{NX}(jac.func(p))
end

function gradient(jac::Union{Linear, Nonlinear}, ::Val{NX}, args...) where NX
    return SVector{NX}(jac.func(args...))
end

# out-of-place, dynamically sized

function gradient(::Empty, nx, ny, args...)
    return Matrix{Float64}(val(nx), val(ny))
end

function gradient(::Zeros, nx, ny, args...)
    return zeros(val(nx), val(ny))
end

function gradient(::Identity, nx, ny, args...)
    return I(min(val(nx), val(ny)))
end

function gradient(grad::Invariant, nx, ny, args...)
    return grad.value
end

function gradient(grad::Constant, nx, ny, p)
    return grad.func(p)
end

function gradient(grad::Constant, nx, ny, dx, x, p, t)
    return grad.func(p)
end

function gradient(grad::Constant, nx, ny, dx, x, y, p, t)
    return grad.func(p)
end

function gradient(grad::Union{Linear, Nonlinear}, nx, ny, args...)
    return grad.func(args...)
end

# in-place

function gradient!(J, ::Empty, args...)
    return J
end

function gradient!(J, ::Zeros, args...)
    return J .= 0
end

function gradient!(J, ::Identity, args...)
    J .= 0
    for i = 1:size(J,1)
        J[i] = 1
    end
    return J
end

function gradient!(J, grad::Invariant, args...)
    return J .= grad.grad
end

function gradient!(J, grad::Constant, p)
    return grad.func(J, p)
end

function gradient!(J, grad::Constant, dx, x, p, t)
    return grad.func(J, p)
end

function gradient!(J, grad::Constant, dx, x, y, p, t)
    return grad.func(J, p)
end

function gradient!(J, grad::Union{Linear, Nonlinear}, args...)
    return grad.func(J, args...)
end


text:
# See https://en.wikipedia.org/wiki/Graph_(discrete_mathematics)

using Julog

clauses = @julog [
    oneway(paris, orleans) <<= true,
    oneway(paris, chartres) <<= true,
    oneway(paris, amiens) <<= true,
    oneway(orleans, blois) <<= true,
    oneway(orleans, bourges) <<= true,
    oneway(blois, tours) <<= true,
    oneway(chartres, lemans) <<= true,
    oneway(lemans, angers) <<= true,
    oneway(lemans, tours) <<= true,
    oneway(angers, nantes) <<= true,

    path(X, Y) <<= oneway(X, Y),
    path(X, Y) <<= oneway(X, Z) & path(Z, Y)
]

cases = [
    :(resolve(@julog(path(X, angers)), clauses)),
    :(resolve(@julog(path(angers, X)), clauses)),
    :(resolve(@julog(path(paris, nantes)), clauses))
]

for case in cases
    println("[] :julia-answer \"\"\"", case, " = ", eval(case), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
"""
    coef(MLMNet)

Extract all coefficients from Mlmnet object

# Arguments

- MLMNet = Mlmnet object

# Value

3d array of coefficients

"""
function coef(MLMNet::Mlmnet)

    return MLMNet.B
end

"""
    coef(MLMNet, lambda)

Extract coefficients from Mlmnet object at a given lambda

# Arguments

- MLMNet = Mlmnet object
- lambda = lambda penalty to use, a floating scalar

# Value

2d array of coefficients

"""
function coef(MLMNet::Mlmnet, lambda::Float64)

    # Find the index corresponding to the lambda of interest
    idx = findall([isapprox(lam, lambda) for lam in MLMNet.lambdas])
    if (length(idx) == 0)
        error("This lambda was not used in the MLMNet object.")
    end

    # Return the coefficient slice for the lambda of interest
    return MLMNet.B[idx[1],:,:]
end

"""
    coef_2d(MLMNet)

Extract coefficients from Mlmnet object as a flattened 2d array

# Arguments

- MLMNet = Mlmnet object

# Value

2d array of flattened coefficients, where each column corresponds to a
different lambda

"""
function coef_2d(MLMNet::Mlmnet)

    # Extract coefficients
    coeffs = coef(MLMNet)

    # Initialize 2d array to store flattened coefficients
    coeffs2d = Array{Float64}(undef, size(coeffs,2)*size(coeffs,3),
                              size(coeffs,1))
    # Iterate through each coefficient slice
    for i in 1:size(coeffs,1)
        coeffs2d[:,i] = vec(coeffs[i,:,:])
    end

    return coeffs2d
end

"""
    predict(MLMNet, lambda, newPredictors)

Calculate new predictions based on Mlmnet object and given a lambda

# Arguments

- MLMNet = Mlmnet object
- lambda = lambda penalty to use, a floating scalar
- newPredictors = Predictors object. Defaults to the data.predictors field
  in the MLM object used to fit the model.

# Value

2d array of predicted values

"""
function predict(MLMNet::Mlmnet, lambda::Float64,
                 newPredictors::Predictors=MLMNet.data.predictors)

    # Include X and Z intercepts in new predictors if necessary
    if MLMNet.data.predictors.isXIntercept==true &&
       newPredictors.isXIntercept==false
        newPredictors.X = add_intercept(newPredictors.X)
        newPredictors.isXIntercept = true
        println("Adding X intercept to newPredictors.")
  	end
    if MLMNet.data.predictors.isZIntercept==true &&
       newPredictors.isZIntercept==false
        newPredictors.Z = add_intercept(newPredictors.Z)
        newPredictors.isZIntercept = true
        println("Adding Z intercept to newPredictors.")
  	end

    # Remove X and Z intercepts in new predictors if necessary
    if MLMNet.data.predictors.isXIntercept==false &&
       newPredictors.isXIntercept==true
        newPredictors.X = remove_intercept(newPredictors.X)
        newPredictors.isXIntercept = false
        println("Removing X intercept from newPredictors.")
    end
    if MLMNet.data.predictors.isZIntercept==false &&
       newPredictors.isZIntercept==true
        newPredictors.Z = remove_intercept(newPredictors.Z)
        newPredictors.isZIntercept = false
        println("Removing Z intercept from newPredictors.")
    end

    # Calculate new predictions
    return calc_preds(newPredictors.X, newPredictors.Z, coef(MLMNet, lambda))
end

"""
    predict(MLMNet, newPredictors)

Calculate new predictions based on Mlmnet object

# Arguments

- MLMNet = Mlmnet object
- newPredictors = Predictors object. Defaults to the data.predictors field
  in the MLM object used to fit the model.

# Value

3d array of predicted values

"""
function predict(MLMNet::Mlmnet,
                 newPredictors::Predictors=MLMNet.data.predictors)

    # Initialize 3d array for storing predictions
    all_preds = Array{Float64}(undef, length(MLMNet.lambdas),
                               size(newPredictors.X,1),
                               size(newPredictors.Z,1))
	# Calculate new predictions for each lambda
    for i = 1:length(MLMNet.lambdas)
        all_preds[i,:,:] = predict(MLMNet, MLMNet.lambdas[i], newPredictors)
    end

    return all_preds
end

"""
    fitted(MLMNet, lambda)

Calculate fitted values of an Mlmnet object, given a lambda

# Arguments

- MLMNet = Mlmnet object
- lambda = lambda penalty to use, a floating scalar

# Value

2d array of fitted values

"""
function fitted(MLMNet::Mlmnet, lambda::Float64)

    return predict(MLMNet, lambda)
end

"""
    fitted(MLMNet)

Calculate fitted values of an Mlmnet object

# Arguments

- MLMNet = Mlmnet object

# Value

3d array of fitted values

"""
function fitted(MLMNet::Mlmnet)

    return predict(MLMNet)
end

"""
    resid(MLMNet, lambda, newData)

Calculate residuals of an MLMNet object, given a lambda

# Arguments

- MLM = Mlmnet object
- lambda = lambda penalty to use, a floating scalar
- newData = RawData object. Defaults to the data field in the MLM object
  used to fit the model.

# Value

2d array of residuals

"""
function resid(MLMNet::Mlmnet, lambda::Float64, newData::RawData=MLMNet.data)

    # Include X and Z intercepts in new data if necessary
    if MLMNet.data.predictors.isXIntercept==true &&
       newData.predictors.isXIntercept==false
        newData.predictors.X = add_intercept(newData.predictors.X)
        newData.predictors.isXIntercept = true
        newData.p = newData.p + 1
        println("Adding X intercept to newData.")
    end
    if MLMNet.data.predictors.isZIntercept==true &&
       newData.predictors.isZIntercept==false
        newData.predictors.Z = add_intercept(newData.predictors.Z)
        newData.predictors.isZIntercept = true
        newData.q = newData.q + 1
        println("Adding Z intercept to newData.")
    end

    # Remove X and Z intercepts in new data if necessary
    if MLMNet.data.predictors.isXIntercept==false &&
       newData.predictors.isXIntercept==true
        newData.predictors.X = remove_intercept(newData.predictors.X)
        newData.predictors.isXIntercept = false
        newData.p = newData.p - 1
        println("Removing X intercept from newData.")
    end
    if MLMNet.data.predictors.isZIntercept==false &&
       newData.predictors.isZIntercept==true
        newData.predictors.Z = remove_intercept(newData.predictors.Z)
        newData.predictors.isZIntercept = false
        newData.q = newData.q - 1
        println("Removing Z intercept from newData.")
    end

    # Calculate residuals
    return calc_resid(get_X(newData), get_Y(newData), get_Z(newData),
                      coef(MLMNet, lambda))
end

"""
    resid(MLMNet, newData)

Calculate residuals of an MLMNet object

# Arguments

- MLMNet = Mlmnet object
- newData = RawData object. Defaults to the data field in the MLM object
  used to fit the model.

# Value

3d array of residuals

"""
function resid(MLMNet::Mlmnet, newData::RawData=MLMNet.data)

    # Initialize 3d array for storing residuals
    all_resid = Array{Float64}(undef, length(MLMNet.lambdas),
                               newData.n, newData.m)

    # Calculate residuals for each lambda
    for i = 1:length(MLMNet.lambdas)
        all_resid[i,:,:] = resid(MLMNet, MLMNet.lambdas[i], newData)
    end

    return all_resid
end


text:
using PMDPs
using PMDPs.LP
using DrWatson

using POMDPs
using BSON, CSV
using DrWatson
using RandomNumbers.Xorshifts
using Random
using DataFrames
using StaticArrays, Distributions # load

using MCTS, DiscreteValueIteration
using POMDPSimulators # load histories

using LightGraphs
using GraphPlot
using Cairo, Compose

using Debugger

using TableView

include(srcdir("MDPPricing.jl"))

objective = :revenue
# pp_params = Dict(pairs((nᵣ=3, c=3, T=10, expected_res=3., res_budget_μ=5., objective=objective)))
pp_params =  Dict(pairs((nᵣ=5, c=3, T=30, expected_res=30., res_budget_μ=5.,
                    objective=objective))) # VI runs
# pp_params =  Dict(pairs((nᵣ=6, c=3, T=30, expected_res=36., res_budget_μ=5.,
#                    objective=objective))) # VI does not run
name = "linear_problem"

"""
Prepare traces
"""

N_individ=10
display("Generating $name with $pp_params")

pp = PMDPs.linear_pp(;pp_params...)
mg = PMDPs.PMDPg(pp)

rnd = Xorshift128Plus(1)
traces = [PMDPs.simulate_trace(mg, rnd) for i in 1:N_individ]

sname = savename("traces_lp", pp_params,  "bson")
@tagsave(datadir("test_traces", sname),  @dict(name, pp, pp_params, traces))

# jldopen(datadir("test_traces", sname), "w") do file
#     addrequire(file, PMDPs)
#     write(file, "data", @strdict(name, pp_params))
# end

# PMDPs.State{3}

# @save(datadir("test_traces", sname),  @strdict(name, pp, pp_params, traces))
# load(datadir("test_traces", sname))

"""
Evaluate
"""
N_sim = 10

dpw_solver = DPWSolver(;depth=50,
                exploration_constant=40.0, max_time=1.,
                enable_state_pw = false,
                keep_tree=true, show_progress=false, rng=Xorshift128Plus())

mcts_solver = MCTSSolver(;depth=50,
                exploration_constant=40.0, max_time=1.,
                rng=Xorshift128Plus())

# pp_params = Dict(pairs((nᵣ=3, c=3, T=10, expected_res=3., res_budget_μ=5., objective=objective)))
# name = "linear_problem"
out_folder="test"

display("Evaluating $name with $pp_params")
sname = savename("traces_lp", pp_params,  "bson")
# data = load(datadir("traces", sname))
# data = load(datadir("test_traces", sname))
data = PMDPs.load_traces(datadir("test_traces", sname))

# PMDPs.process_data(data, PMDPs.flatrate; folder=out_folder, N=N_sim)
# PMDPs.process_data(data, PMDPs.hindsight; folder=out_folder, N=N_sim)
PMDPs.process_data(data, PMDPs.vi; folder=out_folder, N=N_sim)
# PMDPs.process_data(data, PMDPs.fhvi; folder=out_folder, N=N_sim)

# PMDPs.process_data(data, PMDPs.mcts; folder=out_folder, N=N_sim, method_info="dpw", mcts_solver=dpw_solver)
# PMDPs.process_data(data, PMDPs.mcts; folder=out_folder, N=N_sim, method_info="vanilla", mcts_solver=mcts_solver)
println("LP Done.")

"""
Collect results
"""

# results = folder_report(datadir("results", "test", "linear_problem"))
# agg_res = format_result_table(results)

# using WebIO
# vscodedisplay(agg_res)