text:
module WAHTest

using Test
using WAHVectors

function make_random_bitvector(n0, nF, nX, nshuffle)
    v = Vector{UInt32}(n0 + nF + nX)
    aR = rand(0x00000000:0xFFFFFFFF, nX)
    i = 1
    for _ in 1:n0
        v[i] = 0x00000000
        i += 1
    end
    for _ in 1:nF
        v[i] = 0xFFFFFFFF
        i += 1
    end
    for j in 1:nX
        v[i] = aR[j]
    end
    for _ in 1:nshuffle
        shuffle!(v)
    end
    return v
end

@testset "Internals" begin
    @testset "WAH Vector Elements" begin
        @testset "Constructors" begin
            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                @test WAHVectors.WAHElement(0x00000000, n) == convert(WAHVectors.WAHElement, 0x80000000 + n)
                @test WAHVectors.WAHElement(0x00000001, n) == convert(WAHVectors.WAHElement, 0xC0000000 + n)
            end
        end
        @testset "isruns" begin
            @test !WAHVectors.isruns(WAHVectors.WAH_LITERAL_ZEROS)
            @test !WAHVectors.isruns(WAHVectors.WAH_LITERAL_ONES)
            @test WAHVectors.isruns(WAHVectors.WAH_FULL_ZEROS)
            @test WAHVectors.isruns(WAHVectors.WAH_FULL_ONES)

            @test !WAHVectors.is_zeros_runs(WAHVectors.WAH_LITERAL_ZEROS)
            @test !WAHVectors.is_zeros_runs(WAHVectors.WAH_LITERAL_ONES)
            @test WAHVectors.is_zeros_runs(WAHVectors.WAH_FULL_ZEROS)
            @test !WAHVectors.is_zeros_runs(WAHVectors.WAH_FULL_ONES)

            @test !WAHVectors.is_ones_runs(WAHVectors.WAH_LITERAL_ZEROS)
            @test !WAHVectors.is_ones_runs(WAHVectors.WAH_LITERAL_ONES)
            @test !WAHVectors.is_ones_runs(WAHVectors.WAH_FULL_ZEROS)
            @test WAHVectors.is_ones_runs(WAHVectors.WAH_FULL_ONES)

            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e0 = WAHVectors.WAHElement(0x00000000, n)
                e1 = WAHVectors.WAHElement(0x00000001, n)
                el = WAHVectors.WAHElement(rand(0x00000000:0x7FFFFFFF))

                @test WAHVectors.isruns(e0)
                @test WAHVectors.isruns(e1)
                @test !WAHVectors.isruns(el)
                @test WAHVectors.is_zeros_runs(e0)
                @test !WAHVectors.is_zeros_runs(e1)
                @test !WAHVectors.is_zeros_runs(el)
                @test !WAHVectors.is_ones_runs(e0)
                @test WAHVectors.is_ones_runs(e1)
                @test !WAHVectors.is_ones_runs(el)
            end
        end
        @testset "isliteral" begin
            @test WAHVectors.isliteral(WAHVectors.WAH_LITERAL_ZEROS)
            @test WAHVectors.isliteral(WAHVectors.WAH_LITERAL_ONES)
            @test !WAHVectors.isliteral(WAHVectors.WAH_FULL_ZEROS)
            @test !WAHVectors.isliteral(WAHVectors.WAH_FULL_ONES)

            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e0 = WAHVectors.WAHElement(0x00000000, n)
                e1 = WAHVectors.WAHElement(0x00000001, n)
                el = WAHVectors.WAHElement(rand(0x00000000:0x7FFFFFFF))

                @test !WAHVectors.isliteral(e0)
                @test !WAHVectors.isliteral(e1)
                @test WAHVectors.isliteral(el)
            end
        end
        @testset "nwords" begin
            @test WAHVectors.nwords(WAHVectors.WAH_LITERAL_ZEROS) == 1
            @test WAHVectors.nwords(WAHVectors.WAH_LITERAL_ONES) == 1
            @test WAHVectors.nwords(WAHVectors.WAH_FULL_ZEROS) == WAHVectors.WAH_MAX_NRUNS
            @test WAHVectors.nwords(WAHVectors.WAH_FULL_ONES) == WAHVectors.WAH_MAX_NRUNS

            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e0 = WAHVectors.WAHElement(0x00000000, n)
                e1 = WAHVectors.WAHElement(0x00000001, n)
                el = WAHVectors.WAHElement(rand(0x00000000:0x7FFFFFFF))

                @test WAHVectors.nwords(e0) == n
                @test WAHVectors.nwords(e1) == n
                @test WAHVectors.nwords(el) == 0x00000001
            end
        end
        @testset "runval" begin
            @test WAHVectors.runval(WAHVectors.WAH_FULL_ZEROS) == 0x00000000
            @test WAHVectors.runval(WAHVectors.WAH_FULL_ONES) == 0x00000001

            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e0 = WAHVectors.WAHElement(0x00000000, n)
                e1 = WAHVectors.WAHElement(0x00000001, n)

                @test WAHVectors.runval(e0) == 0x00000000
                @test WAHVectors.runval(e1) == 0x00000001
            end
        end
        @testset "nruns" begin
            @test WAHVectors.nruns(WAHVectors.WAH_FULL_ZEROS) == WAHVectors.WAH_MAX_NRUNS
            @test WAHVectors.nruns(WAHVectors.WAH_FULL_ONES) == WAHVectors.WAH_MAX_NRUNS

            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e0 = WAHVectors.WAHElement(0x00000000, n)
                e1 = WAHVectors.WAHElement(0x00000001, n)

                @test WAHVectors.nruns(e0) == n
                @test WAHVectors.nruns(e1) == n
            end
        end
        @testset "nfree" begin
            @test WAHVectors.nfree(WAHVectors.WAH_FULL_ZEROS) == 0
            @test WAHVectors.nfree(WAHVectors.WAH_FULL_ONES) == 0

            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e0 = WAHVectors.WAHElement(0x00000000, n)
                e1 = WAHVectors.WAHElement(0x00000001, n)

                @test WAHVectors.nfree(e0) == (WAHVectors.WAH_MAX_NRUNS - n)
                @test WAHVectors.nfree(e1) == (WAHVectors.WAH_MAX_NRUNS - n)
            end
        end
        @testset "isfull" begin
            @test WAHVectors.isfull(WAHVectors.WAH_FULL_ZEROS)
            @test WAHVectors.isfull(WAHVectors.WAH_FULL_ONES)

            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e0 = WAHVectors.WAHElement(0x00000000, n)
                e1 = WAHVectors.WAHElement(0x00000001, n)

                @test WAHVectors.isfull(e0) == (WAHVectors.WAH_MAX_NRUNS == n)
                @test WAHVectors.isfull(e1) == (WAHVectors.WAH_MAX_NRUNS == n)
            end
        end
        @testset "hasroom" begin
            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                s = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                e = WAHVectors.WAHElement(0x00000000, n)
                @test WAHVectors.hasroom(e) == (n < WAHVectors.WAH_MAX_NRUNS)
                @test WAHVectors.hasroom(e, s) == (s <= WAHVectors.nfree(e))
            end
        end
        @testset "matchingfills" begin
            for _ in 1:100
                n1 = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                n2 = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                v1 = rand(0x00000000:0x00000001)
                v2 = rand(0x00000000:0x00000001)

                @test WAHVectors.matchingfills(WAHVectors.WAHElement(v1, n1), WAHVectors.WAHElement(v2, n2)) == (v1 == v2)
            end
        end
        @testset "increment_nruns_unsafe" begin
            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                i = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                v = rand(0x00000000:0x00000001)
                e = WAHVectors.WAHElement(v, n)
                @test WAHVectors.increment_nruns_unsafe(e) == WAHVectors.WAHElement(UInt32(e) + UInt32(1))
                @test WAHVectors.increment_nruns_unsafe(e, i) == WAHVectors.WAHElement(UInt32(e) + i)
            end
        end
        @testset "runfill" begin
            for _ in 1:100
                n = rand(0x00000001:WAHVectors.WAH_MAX_NRUNS)
                v = rand(0x00000000:0x00000001)
                e1 = WAHVectors.WAHElement(v, n)
                @test WAHVectors.runfill(e1) == ifelse(v == 0x00000001, 0x7FFFFFFF, 0x00000000)
            end
        end
    end
#=
    @testset "WAHVectors" begin
        @testset "round trip" begin
            bv = make_random_bitvector(10, 20, 10, 5)
            @test convert(Vector{UInt32}, WAHVector(bv)) == bv
        end
    end
=#
end

end


text:
#using Pkg
#pwd()
#Pkg.activate(".")
using ChemometricsTools
#View the data in the package space
ChemometricsToolsDatasets()
#Load the data from package-space
DF = ChemometricsToolsDataset("tecator.csv");

Xindices = findall(occursin.( "absorbance", string.(names(DF) ) ) );
names(DF)[(end - 2) : end]
Moisture = DF[ 123 ];
Tecator = convert(Matrix,  DF[:,Xindices]);
#Let's check the Tecator data for duplicate entries ;)
length(AssessHealth(Tecator)["Duplicates"])
#Eeek plenty of duplicates, oh well this is a shoot out not a real analysis...
Df = 0; #Clear the memory...

using Plots
#plot(Tecator', legend = false, ylabel = "Absorbance", xlabel = "bins")

#Shuffle the set, because the classes are ordered
(Tecator, Moisture) = Shuffle(Tecator, Moisture);
#Split the data into a training and hold out set
((TrainX, TrainY), (TestX, TestY)) = SplitByProportion(Tecator, Moisture, 0.7);
msc = MultiplicativeScatterCorrection(TrainX);
TrainX = msc(TrainX);
TestX = msc(TestX);

#Can't have spectra and not plot it!
#plot(TrainX', legend = false, ylabel = "Absorbance", xlabel = "bins")

#Let the shoot out begin!
@time cls   = ClassicLeastSquares( TrainX[:, [5, 25, 45, 60, 90] ], TrainY; Bias = true )
@time rr    = RidgeRegression( TrainX, TrainY, 0.1; Bias = true )
@time krr   = KernelRidgeRegression( TrainX, TrainY, 0.01; KernelParameter = 1.0, KernelType = "gaussian" )
@time lssvm = LSSVM( TrainX, TrainY, 0.01; KernelParameter = 1.0, KernelType = "gaussian" )
pca = PCA(TrainX;Factors = 4);
@time pcr   = PrincipalComponentRegression(pca, TrainY )
@time plsr  = PartialLeastSquares( TrainX, TrainY; Factors = 4, tolerance = 1e-8, maxiters = 200 )
@time elm   = ExtremeLearningMachine(TrainX, TrainY, 11; ActivationFn = sigmoid)
@time rt    = RegressionTree( TrainX, TrainY; gainfn = ssd, maxdepth = 15, minbranchsize = 5)
@time rf    = RandomForest(TrainX, TrainY, :regression; gainfn = ssd, trees = 100,
                       maxdepth = 15,  minbranchsize = 5,
                       samples = 0.7, maxvars = nothing)

#Let's evaluate them on the training set then the hold out - see how they do!
models = ( rr, krr, lssvm, pcr, plsr, elm, rt, rf );
modelnames = ( :rr, :krr, :lssvm, :pcr, :plsr, :elm, :rt, :rf );

TrainingRMSE = Dict()
TrainingRMSE[:cls] = RMSE( cls(TrainX[:,[5, 25, 45, 60, 90]]), TrainY )
for ( name, model ) in zip(modelnames, models)
    TrainingRMSE[name] = RMSE( model(TrainX), TrainY )
    if name == :rf
        println(model(TrainX))
    end
end

TrainingRMSE

HoldOutRMSE = Dict()
HoldOutRMSE[:cls] = RMSE( cls(TestX[:,[5, 25, 45, 60, 90]]), TestY )
for ( name, model ) in zip(modelnames, models)
    HoldOutRMSE[name] = RMSE( model(TestX), TestY )
end
HoldOutRMSE

ext: jl


text:
using Statistics, Random, Distributions

# returns the variable (or matrix), lagged p times,
# with the first p rows filled with ones (to avoid divide errors)
# remember to drop those rows before doing analysis
function lag(x,p::Int64)
	n = size(x,1)
	lagged_x = zeros(eltype(x),n,p)
	lagged_x = [ones(p); x[1:n-p]]
end

# lags of a vector from 1 to p in a p column array
function  lags(x,p)
	n = size(x,1)
	lagged_x = zeros(eltype(x),n,p)
	for i = 1:p
		lagged_x[:,i] = lag(x,i)
	end
    return lagged_x
end

# compute moving average using p most recent values, including current value
function ma(x, p)
    m = similar(x)
    for i = p:size(x,1)
        m[i] = mean(x[i-p+1:i])
    end
    return m
end

# auxiliary model: HAR-RV
# Corsi, Fulvio. "A simple approximate long-memory model
# of realized volatility." Journal of Financial Econometrics 7,
# no. 2 (2009): 174-196.
function HAR(y)
    ylags = lags(y,10)
    X = [ones(size(y,1)) ylags[:,1]  mean(ylags[:,1:4],dims=2) mean(ylags[:,1:10],dims=2)]
    # drop missings
    y = y[11:end]
    X = X[11:end,:]
    βhat = X\y
    σhat = std(y-X*βhat)
    vcat(βhat,σhat)
end

function aux_stat(y)
    y, m, s  = stnorm(abs.(y))
    # look for evidence of volatility clusters
    mm = ma(y,5)
    mm = mm[5:end]
    clusters = quantile(mm,0.75) -quantile(mm, 0.25)
    ϕ = HAR(y)
    vcat(m, s, clusters, ϕ)
end

# the dgp: simple discrete time stochastic volatility (SV) model
function SVmodel(σe, ρ, σu, n, shocks_u, shocks_e)
    burnin = size(shocks_u,1) - n
    hlag = 0.0
    h = ρ.*hlag .+ σu.*shocks_u[1] # figure out type
    y = σe.*exp(h./2.0).*shocks_e[1]
    ys = zeros(eltype(y),n)
    for t = 1:burnin+n
        h = ρ.*hlag .+ σu.*shocks_u[t]
        y = σe.*exp(h./2.0).*shocks_e[t]
        if t > burnin
            ys[t-burnin] = y
        end
        hlag = h
    end
    #plot(ys)
    sqrt(n)*aux_stat(ys)
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule zimg_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("zimg")
JLLWrappers.@generate_main_file("zimg", UUID("1d1c96c1-3ba7-528d-adec-26a13d5085b5"))
end  # module zimg_jll


text:
using POMDPs
using DESPOT
using POMDPToolbox
using POMDPModels

import DESPOT: bounds

struct BabyBounds end
bounds{S}(::BabyBounds, p::BabyPOMDP, ::Vector{DESPOTParticle{S}}, ::DESPOTConfig) = (p.r_feed+p.r_hungry)/(1.0-discount(p)), 0

solver = DESPOTSolver{Bool, Bool, Bool, BabyBounds, MersenneStreamArray}(bounds = BabyBounds(),
                                                    random_streams=MersenneStreamArray(MersenneTwister(1)),
                                                    next_state=false,
                                                    curr_obs=false,
                                                    rng=Base.GLOBAL_RNG
                                                    )
problem = BabyPOMDP()

# test_solver is from POMDPToolbox
test_solver(solver, problem, updater=updater(problem))

test_solver(solver, problem)

struct LightDarkBounds end
bounds{S}(::LightDarkBounds, p::LightDark1D, ::Vector{DESPOTParticle{S}}, ::DESPOTConfig) = p.incorrect_r/(1.0-discount(p)), p.correct_r/(1.0-discount(p))

solver = DESPOTSolver{LightDark1DState, Int64, Float64, LightDarkBounds, MersenneStreamArray}(bounds=LightDarkBounds(),
                                                                  random_streams=MersenneStreamArray(MersenneTwister(1)),
                                                                  next_state=LightDark1DState(0, 0.0),
                                                                  curr_obs=0.0
                                                                 )

test_solver(solver, LightDark1D())


text:
using Test
using Gadgetron
import Base.Threads.@spawn
using Random

testheader = MRD.MRDHeader("""<?xml version="1.0"?>
<ismrmrdHeader xmlns="http://www.ismrm.org/ISMRMRD" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xs="http://www.w3.org/2001/XMLSchema" xsi:schemaLocation="http://www.ismrm.org/ISMRMRD ismrmrd.xsd">
  <subjectInformation>
    <patientName>phantom</patientName>
    <patientWeight_kg>70.3068</patientWeight_kg>
  </subjectInformation>
  <acquisitionSystemInformation>
    <systemVendor>SIEMENS</systemVendor>
    <systemModel>Avanto</systemModel>
    <systemFieldStrength_T>1.494</systemFieldStrength_T>
    <receiverChannels>32</receiverChannels>
    <relativeReceiverNoiseBandwidth>0.79</relativeReceiverNoiseBandwidth>
  </acquisitionSystemInformation>
  <experimentalConditions>
    <H1resonanceFrequency_Hz>63642459</H1resonanceFrequency_Hz>
  </experimentalConditions>
  <encoding>
    <trajectory>cartesian</trajectory>
    <encodedSpace>
      <matrixSize>
        <x>256</x>
        <y>140</y>
        <z>80</z>
      </matrixSize>
      <fieldOfView_mm>
        <x>600</x>
        <y>328.153125</y>
        <z>160</z>
      </fieldOfView_mm>
    </encodedSpace>
    <reconSpace>
      <matrixSize>
        <x>128</x>
        <y>116</y>
        <z>64</z>
      </matrixSize>
      <fieldOfView_mm>
        <x>300</x>
        <y>271.875</y>
        <z>128</z>
      </fieldOfView_mm>
    </reconSpace>
    <encodingLimits>
      <kspace_encoding_step_1>
        <minimum>0</minimum>
        <maximum>83</maximum>
        28
      </kspace_encoding_step_1>
      <kspace_encoding_step_2>
        <minimum>0</minimum>
        <maximum>45</maximum>
        20
      </kspace_encoding_step_2>
      <slice>
        <minimum>0</minimum>
        <maximum>0</maximum>
        0
      </slice>
      <set>
        <minimum>0</minimum>
        <maximum>0</maximum>
        0
      </set>
    </encodingLimits>
    <parallelImaging>
    <accelerationFactor>
      <kspace_encoding_step_1>1</kspace_encoding_step_1>
      <kspace_encoding_step_2>1</kspace_encoding_step_2>
    </accelerationFactor>
    <calibrationMode>other</calibrationMode>
  </parallelImaging>
  </encoding>

  <sequenceParameters>

4.6

    <TE>2.35</TE>
    <TI>300</TI>
  </sequenceParameters>
</ismrmrdHeader>""")

function parrotserver(port)
	t = @spawn begin
		connection = listen(port)
		for msg in connection
			put!(connection,msg)
		end
		close(connection)
	end
	return t

end

function send_and_receive(connection,data )
		put!(connection,data)
		return take!(connection)
end

random_waveform_header() = MRD.WaveformHeader( [rand(Random.TaskLocalRNG(),ftype) for ftype in fieldtypes(MRD.WaveformHeader)]...)

@testset "Connection" begin
	port = 9002
	server = parrotserver(port)
  config = "< _ />"
  header = testheader
  connection = connect("localhost",port,config,header )

  data = randn(Random.TaskLocalRNG(),ComplexF32,(192,32))
	acq = MRD.Acquisition(MRD.AcquisitionHeader(),data)
  @test acq == send_and_receive(connection,acq)

  wav = MRD.Waveform(random_waveform_header(),rand(Random.TaskLocalRNG(),UInt32,(64,8)))
  @test wav == send_and_receive(connection,wav)

  img = MRD.Image(MRD.ImageHeader(),randn(Random.TaskLocalRNG(),ComplexF32,(256,256,3,1)))
  @test img == send_and_receive(connection,img)

  close(connection)
	wait(server)
end


text:
using FileIO: load, save
import LibSndFile
using Plots
pyplot();

run(`afinfo birdland_original.wav`)

sndorg = load("birdland_original.wav")

run(`afconvert -f mp4f birdland_original.wav -o birdland.mp4`);

run(`afinfo birdland.mp4`);

run(`afconvert -f WAVE -d LEI16 birdland.mp4 birdland_aac.wav`);

run(`afinfo birdland_aac.wav`);

sndaac = load("birdland_aac.wav")

sndmp3 = load("birdland_mp3.wav")

xorg = sndorg.data[:,1];
xaac = sndaac.data[:,1];
fs = sndorg.samplerate;
t = (0 : length(xorg)-1) / fs;

snddif = sndorg - sndaac

sndorg - sndmp3


text:
import AoC2019
intcode_program = AoC2019.read_numbers("input21.txt")

intcode_input_short =
map(x->Int(x), collect(
"""NOT T T
AND A T
AND B T
AND C T
NOT T J
AND D J
WALK
"""
))

result = AoC2019.run_program(intcode_program, input=intcode_input_short)
try
  println(join(map(x->Char(x), result.output)))
catch err
  println(result.output)
end

intcode_input_long =
map(x->Int(x), collect(
"""NOT T T
AND A T
AND B T
AND C T
NOT T J
AND D J
OR E T
OR H T
AND T J
RUN
"""
))

result = AoC2019.run_program(intcode_program, input=intcode_input_long)
try
  println(join(map(x->Char(x), result.output)))
catch err
  println(result.output)
end

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
using KhepriRadiance
using Test

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


text:
# SmoothGrad

"""
    smooth_grad(img, model; method=viz_backprop,top_k=1, mean=0, var=1, n_imgs=50, normalized_n_scaled=true)

Method to visually sharpen the gradient based saliency maps. The key idea is to evaluate the average of gradients of images obtained after adding gaussian noise to the original image(`img`). The results suggest the estimated smoothed gradient, leads to visually more coherent sensitivity maps than the unsmoothed gradient, with the resulting visualization aligning better to the human eye with meaningful features.

# Arguments
 - `img`: Original image of interest for visualization
 - `model`: Model to be used for visualization
 - `method=viz_backprop`: Method to be used while evaluating gradients. Can take methods like `viz_guidedbackprop`, `viz_backprop`, `viz_deconvolution` or any other custom function.
 - `top_k=1`: function evaluates the smoothed gradient for `top_k` predictions just as in case of normal `viz_backprop` function.
 - `mean=0`: Mean of the gaussian noise to be added to the image to generate sample space
 - `var=1`: Variance of the gaussian noise to be added to image to generate sample space
 - `n_imgs=50`: Number of images to be generated and used in sample space while evaluation of `smooth_grad`.
 - `normalized_n_scaled=true`: if true, then the noise is first scaled to 0-255 before adding to the image.

Ref: (https://arxiv.org/pdf/1706.03825.pdf)
"""
function smooth_grad(img, model; method=viz_backprop, top_k=1, mean=0, var=1, n_imgs=50, normalized_n_scaled=true)
    smooth_grad = zeros(size(img)...)
    for i in n_imgs
        noise = (var^0.5).*randn(size(img)...) .+ mean
        if (normalized_n_scaled==false)
            noise *= 255.0
        end
        noise = Float32.(noise)
        noised_img = img + noise
        smooth_grad += method(noised_img, model, top_k = top_k)
    end
    smooth_grad/n_imgs
end

# Gradient X Image

"""
    grad_times_image(img, model; top_k=1)

An extension to the saliency approach. Returns an array of products of the `img` with the gradient evaluated with respect to the provided `img` for the `model`, for `top_k` prediction classes.
"""
function grad_times_image(img, model; top_k=1)
    grads = viz_backprop(img, model, top_k=top_k)
    grads_times_image = []
    for i in 1:length(grads)
        push!(grads_times_image, (grads[i][1].*img, grads[i][2],grads[i][3]))
    end
    grads_times_image
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(F^(c*(a+b*x))*(d+e*x)^m, x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*(-1*b*c*e^-1*(d+e*x)*log(F))^(-1m)*(d+e*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^4, x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^4*log(F)^-1+24*F^(c*(a+b*x))*b^-5*c^-5*e^4*log(F)^-5+-24*F^(c*(a+b*x))*b^-4*c^-4*e^3*log(F)^-4*(d+e*x)+-4*e*F^(c*(a+b*x))*b^-2*c^-2*(d+e*x)^3*log(F)^-2+12*F^(c*(a+b*x))*b^-3*c^-3*e^2*(d+e*x)^2*log(F)^-3)
@test integrate(F^(c*(a+b*x))*(d+e*x)^3, x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^3*log(F)^-1+-6*F^(c*(a+b*x))*b^-4*c^-4*e^3*log(F)^-4+-3*e*F^(c*(a+b*x))*b^-2*c^-2*(d+e*x)^2*log(F)^-2+6*F^(c*(a+b*x))*b^-3*c^-3*e^2*log(F)^-3*(d+e*x))
@test integrate(F^(c*(a+b*x))*(d+e*x)^2, x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^2*log(F)^-1+2*F^(c*(a+b*x))*b^-3*c^-3*e^2*log(F)^-3+-2*e*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2*(d+e*x))
@test integrate(F^(c*(a+b*x))*(d+e*x), x) == :(F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1*(d+e*x)+-1*e*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2)
@test integrate(F^(c*(a+b*x)), x) == :(F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1)
@test integrate(F^(c*(a+b*x))*(d+e*x)^-1, x) == :(F^(c*(a+-1*b*d*e^-1))*e^-1*expintegralei(b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-2, x) == :(-1*F^(c*(a+b*x))*e^-1*(d+e*x)^-1+b*c*F^(c*(a+-1*b*d*e^-1))*e^-2*expintegralei(b*c*e^-1*(d+e*x)*log(F))*log(F))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-3, x) == :(-1//2*F^(c*(a+b*x))*e^-1*(d+e*x)^-2+(1/2)*F^(c*(a+-1*b*d*e^-1))*b^2*c^2*e^-3*log(F)^2*expintegralei(b*c*e^-1*(d+e*x)*log(F))+-1//2*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-1*log(F))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-4, x) == :(-1//3*F^(c*(a+b*x))*e^-1*(d+e*x)^-3+-1//6*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-2*log(F)+-1//6*F^(c*(a+b*x))*b^2*c^2*e^-3*(d+e*x)^-1*log(F)^2+1//6*F^(c*(a+-1*b*d*e^-1))*b^3*c^3*e^-4*log(F)^3*expintegralei(b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-5, x) == :(-1//4*F^(c*(a+b*x))*e^-1*(d+e*x)^-4+-1//12*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-3*log(F)+-1//24*F^(c*(a+b*x))*b^2*c^2*e^-3*(d+e*x)^-2*log(F)^2+-1//24*F^(c*(a+b*x))*b^3*c^3*e^-4*(d+e*x)^-1*log(F)^3+1//24*F^(c*(a+-1*b*d*e^-1))*b^4*c^4*e^-5*log(F)^4*expintegralei(b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d^4+e^4*x^4+4*d*e^3*x^3+4*e*x*d^3+6*d^2*e^2*x^2), x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^4*log(F)^-1+24*F^(c*(a+b*x))*b^-5*c^-5*e^4*log(F)^-5+-24*F^(c*(a+b*x))*b^-4*c^-4*e^3*log(F)^-4*(d+e*x)+-4*e*F^(c*(a+b*x))*b^-2*c^-2*(d+e*x)^3*log(F)^-2+12*F^(c*(a+b*x))*b^-3*c^-3*e^2*(d+e*x)^2*log(F)^-3)
@test integrate(F^(c*(a+b*x))*(d^3+e^3*x^3+3*d*e^2*x^2+3*e*x*d^2), x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^3*log(F)^-1+-6*F^(c*(a+b*x))*b^-4*c^-4*e^3*log(F)^-4+-3*e*F^(c*(a+b*x))*b^-2*c^-2*(d+e*x)^2*log(F)^-2+6*F^(c*(a+b*x))*b^-3*c^-3*e^2*log(F)^-3*(d+e*x))
@test integrate(F^(c*(a+b*x))*(d^2+e^2*x^2+2*d*e*x), x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^2*log(F)^-1+2*F^(c*(a+b*x))*b^-3*c^-3*e^2*log(F)^-3+-2*e*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2*(d+e*x))
@test integrate(F^(c*(a+b*x))*(d^2+e^2*x^2+2*d*e*x)^-1, x) == :(-1*F^(c*(a+b*x))*e^-1*(d+e*x)^-1+b*c*F^(c*(a+-1*b*d*e^-1))*e^-2*expintegralei(b*c*e^-1*(d+e*x)*log(F))*log(F))
@test integrate(F^(c*(a+b*x))*(d^3+e^3*x^3+3*d*e^2*x^2+3*e*x*d^2)^-1, x) == :(-1//2*F^(c*(a+b*x))*e^-1*(d+e*x)^-2+(1/2)*F^(c*(a+-1*b*d*e^-1))*b^2*c^2*e^-3*log(F)^2*expintegralei(b*c*e^-1*(d+e*x)*log(F))+-1//2*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-1*log(F))
@test integrate(F^(c*(a+b*x))*(d^4+e^4*x^4+4*d*e^3*x^3+4*e*x*d^3+6*d^2*e^2*x^2)^-1, x) == :(-1//3*F^(c*(a+b*x))*e^-1*(d+e*x)^-3+-1//6*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-2*log(F)+-1//6*F^(c*(a+b*x))*b^2*c^2*e^-3*(d+e*x)^-1*log(F)^2+1//6*F^(c*(a+-1*b*d*e^-1))*b^3*c^3*e^-4*log(F)^3*expintegralei(b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d^5+e^5*x^5+5*d*e^4*x^4+5*e*x*d^4+10*d^2*e^3*x^3+10*d^3*e^2*x^2)^-1, x) == :(-1//4*F^(c*(a+b*x))*e^-1*(d+e*x)^-4+-1//12*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-3*log(F)+-1//24*F^(c*(a+b*x))*b^2*c^2*e^-3*(d+e*x)^-2*log(F)^2+-1//24*F^(c*(a+b*x))*b^3*c^3*e^-4*(d+e*x)^-1*log(F)^3+1//24*F^(c*(a+-1*b*d*e^-1))*b^4*c^4*e^-5*log(F)^4*expintegralei(b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*((d+e*x)^n)^m, x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*((d+e*x)^n)^m*(-1*b*c*e^-1*(d+e*x)*log(F))^(-1*m*n)*log(F)^-1*SpecialFunctions.gamma(1+m*n,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d^4+e^4*x^4+4*d*e^3*x^3+4*e*x*d^3+6*d^2*e^2*x^2)^m, x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*((d+e*x)^4)^m*(-1*b*c*e^-1*(d+e*x)*log(F))^(-4m)*log(F)^-1*SpecialFunctions.gamma(1+4m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d^3+e^3*x^3+3*d*e^2*x^2+3*e*x*d^2)^m, x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*((d+e*x)^3)^m*(-1*b*c*e^-1*(d+e*x)*log(F))^(-3m)*log(F)^-1*SpecialFunctions.gamma(1+3m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d^2+e^2*x^2+2*d*e*x)^m, x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*((d+e*x)^2)^m*(-1*b*c*e^-1*(d+e*x)*log(F))^(-2m)*log(F)^-1*SpecialFunctions.gamma(1+2m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^m, x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*(-1*b*c*e^-1*(d+e*x)*log(F))^(-1m)*(d+e*x)^m*log(F)^-1*SpecialFunctions.gamma(1+m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^(-1m), x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*(-1*b*c*e^-1*(d+e*x)*log(F))^m*(d+e*x)^(-1m)*log(F)^-1*SpecialFunctions.gamma(1+-1m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d^2+e^2*x^2+2*d*e*x)^(-1m), x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*((d+e*x)^2)^(-1m)*(-1*b*c*e^-1*(d+e*x)*log(F))^(2m)*log(F)^-1*SpecialFunctions.gamma(1+-2m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d^3+e^3*x^3+3*d*e^2*x^2+3*e*x*d^2)^(-1m), x) == :(F^(c*(a+-1*b*d*e^-1))*b^-1*c^-1*((d+e*x)^3)^(-1m)*(-1*b*c*e^-1*(d+e*x)*log(F))^(3m)*log(F)^-1*SpecialFunctions.gamma(1+-3m,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate(F^(2+5x), x) == :(1//5*F^(2+5x)*log(F)^-1)
@test integrate(F^(a+b*x), x) == :(F^(a+b*x)*b^-1*log(F)^-1)
@test integrate(10^(2+5x), x) == :(2^(2+5x)*5^(1+5x)*log(10)^-1)
@test integrate(F^(a+b*x)*x^7//2, x) == :(F^(a+b*x)*b^-1*x^7//2*log(F)^-1+-105//8*F^(a+b*x)*b^-4*x^(1/2)*log(F)^-4+-7//2*F^(a+b*x)*b^-2*x^5//2*log(F)^-2+35//4*F^(a+b*x)*b^-3*x^3//2*log(F)^-3+105//16*pi^(1/2)*F^a*b^-9//2*log(F)^-9//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^5//2, x) == :(F^(a+b*x)*b^-1*x^5//2*log(F)^-1+-5//2*F^(a+b*x)*b^-2*x^3//2*log(F)^-2+15//4*F^(a+b*x)*b^-3*x^(1/2)*log(F)^-3+-15//8*pi^(1/2)*F^a*b^-7//2*log(F)^-7//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^3//2, x) == :(F^(a+b*x)*b^-1*x^3//2*log(F)^-1+-3//2*F^(a+b*x)*b^-2*x^(1/2)*log(F)^-2+3//4*pi^(1/2)*F^a*b^-5//2*log(F)^-5//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^(1/2), x) == :(F^(a+b*x)*b^-1*x^(1/2)*log(F)^-1+-1//2*pi^(1/2)*F^a*b^-3//2*log(F)^-3//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^-1//2, x) == :(pi^(1/2)*F^a*b^-1//2*log(F)^-1//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^-3//2, x) == :(-2*F^(a+b*x)*x^-1//2+2*pi^(1/2)*F^a*b^(1/2)*log(F)^(1/2)*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^-5//2, x) == :(-2//3*F^(a+b*x)*x^-3//2+-4//3*b*F^(a+b*x)*x^-1//2*log(F)+4//3*pi^(1/2)*F^a*b^3//2*log(F)^3//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^-7//2, x) == :(-2//5*F^(a+b*x)*x^-5//2+-8//15*F^(a+b*x)*b^2*x^-1//2*log(F)^2+-4//15*b*F^(a+b*x)*x^-3//2*log(F)+8//15*pi^(1/2)*F^a*b^5//2*log(F)^5//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(a+b*x)*x^-9//2, x) == :(-2//7*F^(a+b*x)*x^-7//2+-16//105*F^(a+b*x)*b^3*x^-1//2*log(F)^3+-8//105*F^(a+b*x)*b^2*x^-3//2*log(F)^2+-4//35*b*F^(a+b*x)*x^-5//2*log(F)+16//105*pi^(1/2)*F^a*b^7//2*log(F)^7//2*erfi(b^(1/2)*x^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^7//2, x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^7//2*log(F)^-1+-105//8*F^(c*(a+b*x))*b^-4*c^-4*e^3*(d+e*x)^(1/2)*log(F)^-4+-7//2*e*F^(c*(a+b*x))*b^-2*c^-2*(d+e*x)^5//2*log(F)^-2+35//4*F^(c*(a+b*x))*b^-3*c^-3*e^2*(d+e*x)^3//2*log(F)^-3+105//16*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^-9//2*c^-9//2*e^7//2*log(F)^-9//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^5//2, x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^5//2*log(F)^-1+-5//2*e*F^(c*(a+b*x))*b^-2*c^-2*(d+e*x)^3//2*log(F)^-2+15//4*F^(c*(a+b*x))*b^-3*c^-3*e^2*(d+e*x)^(1/2)*log(F)^-3+-15//8*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^-7//2*c^-7//2*e^5//2*log(F)^-7//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^3//2, x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^3//2*log(F)^-1+-3//2*e*F^(c*(a+b*x))*b^-2*c^-2*(d+e*x)^(1/2)*log(F)^-2+3//4*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^-5//2*c^-5//2*e^3//2*log(F)^-5//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^(1/2), x) == :(F^(c*(a+b*x))*b^-1*c^-1*(d+e*x)^(1/2)*log(F)^-1+-1//2*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^-3//2*c^-3//2*e^(1/2)*log(F)^-3//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-1//2, x) == :(pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^-1//2*c^-1//2*e^-1//2*log(F)^-1//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-3//2, x) == :(-2*F^(c*(a+b*x))*e^-1*(d+e*x)^-1//2+2*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^(1/2)*c^(1/2)*e^-3//2*log(F)^(1/2)*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-5//2, x) == :(-2//3*F^(c*(a+b*x))*e^-1*(d+e*x)^-3//2+-4//3*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-1//2*log(F)+4//3*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^3//2*c^3//2*e^-5//2*log(F)^3//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-7//2, x) == :(-2//5*F^(c*(a+b*x))*e^-1*(d+e*x)^-5//2+-8//15*F^(c*(a+b*x))*b^2*c^2*e^-3*(d+e*x)^-1//2*log(F)^2+-4//15*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-3//2*log(F)+8//15*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^5//2*c^5//2*e^-7//2*log(F)^5//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^-9//2, x) == :(-2//7*F^(c*(a+b*x))*e^-1*(d+e*x)^-7//2+-16//105*F^(c*(a+b*x))*b^3*c^3*e^-4*(d+e*x)^-1//2*log(F)^3+-8//105*F^(c*(a+b*x))*b^2*c^2*e^-3*(d+e*x)^-3//2*log(F)^2+-4//35*b*c*F^(c*(a+b*x))*e^-2*(d+e*x)^-5//2*log(F)+16//105*pi^(1/2)*F^(c*(a+-1*b*d*e^-1))*b^7//2*c^7//2*e^-9//2*log(F)^7//2*erfi(b^(1/2)*c^(1/2)*e^-1//2*(d+e*x)^(1/2)*log(F)^(1/2)))
@test integrate(x^13//2*exp(-1*b*x), x) == :(-1*b^-1*x^13//2*exp(-1*b*x)+-135135//64*b^-7*x^(1/2)*exp(-1*b*x)+-45045//32*b^-6*x^3//2*exp(-1*b*x)+-9009//16*b^-5*x^5//2*exp(-1*b*x)+-1287//8*b^-4*x^7//2*exp(-1*b*x)+-143//4*b^-3*x^9//2*exp(-1*b*x)+-13//2*b^-2*x^11//2*exp(-1*b*x)+135135//128*pi^(1/2)*b^-15//2*erf(b^(1/2)*x^(1/2)))
@test integrate(F^(c*(a+b*x))*(d+e*x)^4//3, x) == :(-1*e*F^(c*(a+-1*b*d*e^-1))*b^-2*c^-2*(-1*b*c*e^-1*(d+e*x)*log(F))^-1//3*(d+e*x)^1//3*log(F)^-2*SpecialFunctions.gamma(7//3,-1*b*c*e^-1*(d+e*x)*log(F)))
@test integrate((F^(c*(a+b*x)))^n*(d+e*x)^4//3, x) == :(-1*e*F^(c*n*(a+-1*b*d*e^-1)+-1*c*n*(a+b*x))*b^-2*c^-2*n^-2*(F^(c*(a+b*x)))^n*(-1*b*c*n*e^-1*(d+e*x)*log(F))^-1//3*(d+e*x)^1//3*log(F)^-2*SpecialFunctions.gamma(7//3,-1*b*c*n*e^-1*(d+e*x)*log(F)))
@test integrate(F^(c*(a+b*x))*(d+e*x), x) == :(F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1*(d+e*x)+-1*e*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2)
@test integrate(F^(c*(a+b*x))*(d+e*x+f*x^2), x) == :(d*F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1+-1*e*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2+2*f*F^(c*(a+b*x))*b^-3*c^-3*log(F)^-3+e*x*F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1+f*F^(c*(a+b*x))*b^-1*c^-1*x^2*log(F)^-1+-2*f*x*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2)
@test integrate(F^(c*(a+b*x))*(d+e*x+f*x^2+g*x^3), x) == :(d*F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1+-1*e*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2+-6*g*F^(c*(a+b*x))*b^-4*c^-4*log(F)^-4+2*f*F^(c*(a+b*x))*b^-3*c^-3*log(F)^-3+e*x*F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1+f*F^(c*(a+b*x))*b^-1*c^-1*x^2*log(F)^-1+g*F^(c*(a+b*x))*b^-1*c^-1*x^3*log(F)^-1+-3*g*F^(c*(a+b*x))*b^-2*c^-2*x^2*log(F)^-2+-2*f*x*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2+6*g*x*F^(c*(a+b*x))*b^-3*c^-3*log(F)^-3)
@test integrate(F^(c*(a+b*x))*(d+e*x+f*x^2+g*x^3+h*x^4), x) == :(d*F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1+-1*e*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2+-6*g*F^(c*(a+b*x))*b^-4*c^-4*log(F)^-4+2*f*F^(c*(a+b*x))*b^-3*c^-3*log(F)^-3+24*h*F^(c*(a+b*x))*b^-5*c^-5*log(F)^-5+e*x*F^(c*(a+b*x))*b^-1*c^-1*log(F)^-1+f*F^(c*(a+b*x))*b^-1*c^-1*x^2*log(F)^-1+g*F^(c*(a+b*x))*b^-1*c^-1*x^3*log(F)^-1+h*F^(c*(a+b*x))*b^-1*c^-1*x^4*log(F)^-1+-24*h*x*F^(c*(a+b*x))*b^-4*c^-4*log(F)^-4+-4*h*F^(c*(a+b*x))*b^-2*c^-2*x^3*log(F)^-2+-3*g*F^(c*(a+b*x))*b^-2*c^-2*x^2*log(F)^-2+-2*f*x*F^(c*(a+b*x))*b^-2*c^-2*log(F)^-2+6*g*x*F^(c*(a+b*x))*b^-3*c^-3*log(F)^-3+12*h*F^(c*(a+b*x))*b^-3*c^-3*x^2*log(F)^-3)
@test integrate(x^m*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(-1*b^-1*x^m*(b*x)^(-1m)*SpecialFunctions.gamma(4+m,b*x)*exp(-1a)+-1*a^3*b^-1*x^m*(b*x)^(-1m)*SpecialFunctions.gamma(1+m,b*x)*exp(-1a)+-3*a*b^-1*x^m*(b*x)^(-1m)*SpecialFunctions.gamma(3+m,b*x)*exp(-1a)+-3*a^2*b^-1*x^m*(b*x)^(-1m)*SpecialFunctions.gamma(2+m,b*x)*exp(-1a))
@test integrate(x^3*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(-720*b^-4*exp(-1a+-1*b*x)+-30*x^4*exp(-1a+-1*b*x)+-1*b^2*x^6*exp(-1a+-1*b*x)+-720*x*b^-3*exp(-1a+-1*b*x)+-360*a*b^-4*exp(-1a+-1*b*x)+-360*b^-2*x^2*exp(-1a+-1*b*x)+-120*b^-1*x^3*exp(-1a+-1*b*x)+-72*a^2*b^-4*exp(-1a+-1*b*x)+-15*a*x^4*exp(-1a+-1*b*x)+-6*b*x^5*exp(-1a+-1*b*x)+-6*a^3*b^-4*exp(-1a+-1*b*x)+-3*a^2*x^4*exp(-1a+-1*b*x)+-1*a^3*b^-1*x^3*exp(-1a+-1*b*x)+-360*a*x*b^-3*exp(-1a+-1*b*x)+-180*a*b^-2*x^2*exp(-1a+-1*b*x)+-72*x*a^2*b^-3*exp(-1a+-1*b*x)+-60*a*b^-1*x^3*exp(-1a+-1*b*x)+-36*a^2*b^-2*x^2*exp(-1a+-1*b*x)+-12*a^2*b^-1*x^3*exp(-1a+-1*b*x)+-6*x*a^3*b^-3*exp(-1a+-1*b*x)+-3*a*b*x^5*exp(-1a+-1*b*x)+-3*a^3*b^-2*x^2*exp(-1a+-1*b*x))
@test integrate(x^2*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(-120*b^-3*exp(-1a+-1*b*x)+-20*x^3*exp(-1a+-1*b*x)+-1*b^2*x^5*exp(-1a+-1*b*x)+-120*x*b^-2*exp(-1a+-1*b*x)+-72*a*b^-3*exp(-1a+-1*b*x)+-60*b^-1*x^2*exp(-1a+-1*b*x)+-18*a^2*b^-3*exp(-1a+-1*b*x)+-12*a*x^3*exp(-1a+-1*b*x)+-5*b*x^4*exp(-1a+-1*b*x)+-3*a^2*x^3*exp(-1a+-1*b*x)+-2*a^3*b^-3*exp(-1a+-1*b*x)+-1*a^3*b^-1*x^2*exp(-1a+-1*b*x)+-72*a*x*b^-2*exp(-1a+-1*b*x)+-36*a*b^-1*x^2*exp(-1a+-1*b*x)+-18*x*a^2*b^-2*exp(-1a+-1*b*x)+-9*a^2*b^-1*x^2*exp(-1a+-1*b*x)+-3*a*b*x^4*exp(-1a+-1*b*x)+-2*x*a^3*b^-2*exp(-1a+-1*b*x))
@test integrate(x*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(-24*b^-2*exp(-1a+-1*b*x)+-1*b^-2*(a+b*x)^4*exp(-1a+-1*b*x)+-24*b^-2*(a+b*x)*exp(-1a+-1*b*x)+-12*b^-2*(a+b*x)^2*exp(-1a+-1*b*x)+-4*b^-2*(a+b*x)^3*exp(-1a+-1*b*x)+6*a*b^-2*exp(-1a+-1*b*x)+a*b^-2*(a+b*x)^3*exp(-1a+-1*b*x)+3*a*b^-2*(a+b*x)^2*exp(-1a+-1*b*x)+6*a*b^-2*(a+b*x)*exp(-1a+-1*b*x))
@test integrate((a+b*x)^3*exp(-1a+-1*b*x), x) == :(-6*b^-1*exp(-1a+-1*b*x)+-1*b^-1*(a+b*x)^3*exp(-1a+-1*b*x)+-6*b^-1*(a+b*x)*exp(-1a+-1*b*x)+-3*b^-1*(a+b*x)^2*exp(-1a+-1*b*x))
@test integrate(x^-1*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(-2*exp(-1a+-1*b*x)+-3*a*exp(-1a+-1*b*x)+-3*a^2*exp(-1a+-1*b*x)+a^3*exp(-1a)*expintegralei(-1*b*x)+-1*b^2*x^2*exp(-1a+-1*b*x)+-2*b*x*exp(-1a+-1*b*x)+-3*a*b*x*exp(-1a+-1*b*x))
@test integrate(x^-2*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(-1*b*exp(-1a+-1*b*x)+-1*x*b^2*exp(-1a+-1*b*x)+-1*a^3*x^-1*exp(-1a+-1*b*x)+-3*a*b*exp(-1a+-1*b*x)+-1*b*a^3*exp(-1a)*expintegralei(-1*b*x)+3*b*a^2*exp(-1a)*expintegralei(-1*b*x))
@test integrate(x^-3*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(-1*b^2*exp(-1a+-1*b*x)+-1//2*a^3*x^-2*exp(-1a+-1*b*x)+(1/2)*b*a^3*x^-1*exp(-1a+-1*b*x)+(1/2)*a^3*b^2*exp(-1a)*expintegralei(-1*b*x)+-3*b*a^2*x^-1*exp(-1a+-1*b*x)+-3*a^2*b^2*exp(-1a)*expintegralei(-1*b*x)+3*a*b^2*exp(-1a)*expintegralei(-1*b*x))
@test integrate(x^-4*(a+b*x)^3*exp(-1a+-1*b*x), x) == :(b^3*exp(-1a)*expintegralei(-1*b*x)+-1//3*a^3*x^-3*exp(-1a+-1*b*x)+-3*a*b^2*x^-1*exp(-1a+-1*b*x)+-3*a*b^3*exp(-1a)*expintegralei(-1*b*x)+-3//2*b*a^2*x^-2*exp(-1a+-1*b*x)+-1//6*a^3*b^2*x^-1*exp(-1a+-1*b*x)+-1//6*a^3*b^3*exp(-1a)*expintegralei(-1*b*x)+1//6*b*a^3*x^-2*exp(-1a+-1*b*x)+3//2*a^2*b^2*x^-1*exp(-1a+-1*b*x)+3//2*a^2*b^3*exp(-1a)*expintegralei(-1*b*x))
@test integrate(F^(a+b*(c+d*x))*x^m*(e+f*x)^2, x) == :(F^(a+b*c)*b^-1*d^-1*e^2*x^m*(-1*b*d*x*log(F))^(-1m)*log(F)^-1*SpecialFunctions.gamma(1+m,-1*b*d*x*log(F))+F^(a+b*c)*b^-3*d^-3*f^2*x^m*(-1*b*d*x*log(F))^(-1m)*log(F)^-3*SpecialFunctions.gamma(3+m,-1*b*d*x*log(F))+-2*e*f*F^(a+b*c)*b^-2*d^-2*x^m*(-1*b*d*x*log(F))^(-1m)*log(F)^-2*SpecialFunctions.gamma(2+m,-1*b*d*x*log(F)))
@test integrate(F^(a+b*(c+d*x))*x^3*(e+f*x)^2, x) == :(-120*F^(a+b*c+b*d*x)*b^-6*d^-6*f^2*log(F)^-6+-6*F^(a+b*c+b*d*x)*b^-4*d^-4*e^2*log(F)^-4+F^(a+b*c+b*d*x)*b^-1*d^-1*e^2*x^3*log(F)^-1+F^(a+b*c+b*d*x)*b^-1*d^-1*f^2*x^5*log(F)^-1+-60*F^(a+b*c+b*d*x)*b^-4*d^-4*f^2*x^2*log(F)^-4+-5*F^(a+b*c+b*d*x)*b^-2*d^-2*f^2*x^4*log(F)^-2+-3*F^(a+b*c+b*d*x)*b^-2*d^-2*e^2*x^2*log(F)^-2+6*x*F^(a+b*c+b*d*x)*b^-3*d^-3*e^2*log(F)^-3+20*F^(a+b*c+b*d*x)*b^-3*d^-3*f^2*x^3*log(F)^-3+48*e*f*F^(a+b*c+b*d*x)*b^-5*d^-5*log(F)^-5+120*x*F^(a+b*c+b*d*x)*b^-5*d^-5*f^2*log(F)^-5+-48*e*f*x*F^(a+b*c+b*d*x)*b^-4*d^-4*log(F)^-4+-8*e*f*F^(a+b*c+b*d*x)*b^-2*d^-2*x^3*log(F)^-2+2*e*f*F^(a+b*c+b*d*x)*b^-1*d^-1*x^4*log(F)^-1+24*e*f*F^(a+b*c+b*d*x)*b^-3*d^-3*x^2*log(F)^-3)
@test integrate(F^(a+b*(c+d*x))*x^2*(e+f*x)^2, x) == :(2*F^(a+b*c+b*d*x)*b^-3*d^-3*e^2*log(F)^-3+24*F^(a+b*c+b*d*x)*b^-5*d^-5*f^2*log(F)^-5+F^(a+b*c+b*d*x)*b^-1*d^-1*e^2*x^2*log(F)^-1+F^(a+b*c+b*d*x)*b^-1*d^-1*f^2*x^4*log(F)^-1+-24*x*F^(a+b*c+b*d*x)*b^-4*d^-4*f^2*log(F)^-4+-12*e*f*F^(a+b*c+b*d*x)*b^-4*d^-4*log(F)^-4+-4*F^(a+b*c+b*d*x)*b^-2*d^-2*f^2*x^3*log(F)^-2+-2*x*F^(a+b*c+b*d*x)*b^-2*d^-2*e^2*log(F)^-2+12*F^(a+b*c+b*d*x)*b^-3*d^-3*f^2*x^2*log(F)^-3+-6*e*f*F^(a+b*c+b*d*x)*b^-2*d^-2*x^2*log(F)^-2+2*e*f*F^(a+b*c+b*d*x)*b^-1*d^-1*x^3*log(F)^-1+12*e*f*x*F^(a+b*c+b*d*x)*b^-3*d^-3*log(F)^-3)
@test integrate(x*F^(a+b*(c+d*x))*(e+f*x)^2, x) == :(-1*F^(a+b*c+b*d*x)*b^-2*d^-2*e^2*log(F)^-2+-6*F^(a+b*c+b*d*x)*b^-4*d^-4*f^2*log(F)^-4+x*F^(a+b*c+b*d*x)*b^-1*d^-1*e^2*log(F)^-1+F^(a+b*c+b*d*x)*b^-1*d^-1*f^2*x^3*log(F)^-1+-3*F^(a+b*c+b*d*x)*b^-2*d^-2*f^2*x^2*log(F)^-2+4*e*f*F^(a+b*c+b*d*x)*b^-3*d^-3*log(F)^-3+6*x*F^(a+b*c+b*d*x)*b^-3*d^-3*f^2*log(F)^-3+-4*e*f*x*F^(a+b*c+b*d*x)*b^-2*d^-2*log(F)^-2+2*e*f*F^(a+b*c+b*d*x)*b^-1*d^-1*x^2*log(F)^-1)
@test integrate(F^(a+b*(c+d*x))*(e+f*x)^2, x) == :(F^(a+b*c+b*d*x)*b^-1*d^-1*(e+f*x)^2*log(F)^-1+2*F^(a+b*c+b*d*x)*b^-3*d^-3*f^2*log(F)^-3+-2*f*F^(a+b*c+b*d*x)*b^-2*d^-2*log(F)^-2*(e+f*x))
@test integrate(F^(a+b*(c+d*x))*x^-1*(e+f*x)^2, x) == :(F^(a+b*c)*e^2*expintegralei(b*d*x*log(F))+-1*F^(a+b*c+b*d*x)*b^-2*d^-2*f^2*log(F)^-2+x*F^(a+b*c+b*d*x)*b^-1*d^-1*f^2*log(F)^-1+2*e*f*F^(a+b*c+b*d*x)*b^-1*d^-1*log(F)^-1)
@test integrate(F^(a+b*(c+d*x))*x^-2*(e+f*x)^2, x) == :(-1*F^(a+b*c+b*d*x)*e^2*x^-1+2*e*f*F^(a+b*c)*expintegralei(b*d*x*log(F))+F^(a+b*c+b*d*x)*b^-1*d^-1*f^2*log(F)^-1+b*d*F^(a+b*c)*e^2*expintegralei(b*d*x*log(F))*log(F))
@test integrate(F^(a+b*(c+d*x))*x^-3*(e+f*x)^2, x) == :(F^(a+b*c)*f^2*expintegralei(b*d*x*log(F))+-1//2*F^(a+b*c+b*d*x)*e^2*x^-2+-2*e*f*F^(a+b*c+b*d*x)*x^-1+(1/2)*F^(a+b*c)*b^2*d^2*e^2*log(F)^2*expintegralei(b*d*x*log(F))+-1//2*b*d*F^(a+b*c+b*d*x)*e^2*x^-1*log(F)+2*b*d*e*f*F^(a+b*c)*expintegralei(b*d*x*log(F))*log(F))
@test integrate(F^(a+b*(c+d*x))*x^-4*(e+f*x)^2, x) == :(-1*F^(a+b*c+b*d*x)*f^2*x^-1+-1//3*F^(a+b*c+b*d*x)*e^2*x^-3+-1*e*f*F^(a+b*c+b*d*x)*x^-2+b*d*F^(a+b*c)*f^2*expintegralei(b*d*x*log(F))*log(F)+-1//6*b*d*F^(a+b*c+b*d*x)*e^2*x^-2*log(F)+-1//6*F^(a+b*c+b*d*x)*b^2*d^2*e^2*x^-1*log(F)^2+1//6*F^(a+b*c)*b^3*d^3*e^2*log(F)^3*expintegralei(b*d*x*log(F))+e*f*F^(a+b*c)*b^2*d^2*log(F)^2*expintegralei(b*d*x*log(F))+-1*b*d*e*f*F^(a+b*c+b*d*x)*x^-1*log(F))
@test integrate(F^(a+b*(c+d*x))*x^-5*(e+f*x)^2, x) == :(-1//2*F^(a+b*c+b*d*x)*f^2*x^-2+-1//4*F^(a+b*c+b*d*x)*e^2*x^-4+-2//3*e*f*F^(a+b*c+b*d*x)*x^-3+(1/2)*F^(a+b*c)*b^2*d^2*f^2*log(F)^2*expintegralei(b*d*x*log(F))+-1//2*b*d*F^(a+b*c+b*d*x)*f^2*x^-1*log(F)+-1//12*b*d*F^(a+b*c+b*d*x)*e^2*x^-3*log(F)+-1//24*F^(a+b*c+b*d*x)*b^2*d^2*e^2*x^-2*log(F)^2+-1//24*F^(a+b*c+b*d*x)*b^3*d^3*e^2*x^-1*log(F)^3+1//24*F^(a+b*c)*b^4*d^4*e^2*log(F)^4*expintegralei(b*d*x*log(F))+-1//3*b*d*e*f*F^(a+b*c+b*d*x)*x^-2*log(F)+-1//3*e*f*F^(a+b*c+b*d*x)*b^2*d^2*x^-1*log(F)^2+1//3*e*f*F^(a+b*c)*b^3*d^3*log(F)^3*expintegralei(b*d*x*log(F)))
@test integrate((a+b*x)^4*(c+d*x)^3*exp(-1a+-1*b*x), x) == :(-5040*b^-4*d^3*exp(-1a+-1*b*x)+-24*b^-4*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+-1*b^-4*d^3*(a+b*x)^7*exp(-1a+-1*b*x)+-1*b^-4*(a+b*x)^4*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+-5040*b^-4*d^3*(a+b*x)*exp(-1a+-1*b*x)+-2520*b^-4*d^3*(a+b*x)^2*exp(-1a+-1*b*x)+-2160*b^-4*d^2*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-840*b^-4*d^3*(a+b*x)^3*exp(-1a+-1*b*x)+-360*d*b^-4*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-210*b^-4*d^3*(a+b*x)^4*exp(-1a+-1*b*x)+-42*b^-4*d^3*(a+b*x)^5*exp(-1a+-1*b*x)+-24*b^-4*(b*c+-1*a*d)^3*(a+b*x)*exp(-1a+-1*b*x)+-12*b^-4*(a+b*x)^2*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+-7*b^-4*d^3*(a+b*x)^6*exp(-1a+-1*b*x)+-4*b^-4*(a+b*x)^3*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+-2160*b^-4*d^2*(a+b*x)*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-1080*b^-4*d^2*(a+b*x)^2*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-360*d*b^-4*(b*c+-1*a*d)^2*(a+b*x)*exp(-1a+-1*b*x)+-360*b^-4*d^2*(a+b*x)^3*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-180*d*b^-4*(a+b*x)^2*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-90*b^-4*d^2*(a+b*x)^4*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-60*d*b^-4*(a+b*x)^3*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-18*b^-4*d^2*(a+b*x)^5*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-15*d*b^-4*(a+b*x)^4*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-3*d*b^-4*(a+b*x)^5*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-3*b^-4*d^2*(a+b*x)^6*(b*c+-1*a*d)*exp(-1a+-1*b*x))
@test integrate((a+b*x)^4*(c+d*x)^2*exp(-1a+-1*b*x), x) == :(-720*b^-3*d^2*exp(-1a+-1*b*x)+-24*b^-3*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-1*b^-3*d^2*(a+b*x)^6*exp(-1a+-1*b*x)+-1*b^-3*(a+b*x)^4*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-720*b^-3*d^2*(a+b*x)*exp(-1a+-1*b*x)+-360*b^-3*d^2*(a+b*x)^2*exp(-1a+-1*b*x)+-240*d*b^-3*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-120*b^-3*d^2*(a+b*x)^3*exp(-1a+-1*b*x)+-30*b^-3*d^2*(a+b*x)^4*exp(-1a+-1*b*x)+-24*b^-3*(b*c+-1*a*d)^2*(a+b*x)*exp(-1a+-1*b*x)+-12*b^-3*(a+b*x)^2*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-6*b^-3*d^2*(a+b*x)^5*exp(-1a+-1*b*x)+-4*b^-3*(a+b*x)^3*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-240*d*b^-3*(a+b*x)*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-120*d*b^-3*(a+b*x)^2*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-40*d*b^-3*(a+b*x)^3*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-10*d*b^-3*(a+b*x)^4*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-2*d*b^-3*(a+b*x)^5*(b*c+-1*a*d)*exp(-1a+-1*b*x))
@test integrate((a+b*x)^4*(c+d*x)*exp(-1a+-1*b*x), x) == :(-1*b^-2*(-24*a*d+24*b*c)*exp(-1a+-1*b*x)+-120*d*b^-2*exp(-1a+-1*b*x)+-1*d*b^-2*(a+b*x)^5*exp(-1a+-1*b*x)+-1*b^-2*(a+b*x)^2*(-12*a*d+12*b*c)*exp(-1a+-1*b*x)+-1*b^-2*(a+b*x)^3*(-4*a*d+4*b*c)*exp(-1a+-1*b*x)+-1*b^-2*(a+b*x)^4*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-1*b^-2*(a+b*x)*(-24*a*d+24*b*c)*exp(-1a+-1*b*x)+-120*d*b^-2*(a+b*x)*exp(-1a+-1*b*x)+-60*d*b^-2*(a+b*x)^2*exp(-1a+-1*b*x)+-20*d*b^-2*(a+b*x)^3*exp(-1a+-1*b*x)+-5*d*b^-2*(a+b*x)^4*exp(-1a+-1*b*x))
@test integrate((a+b*x)^4*exp(-1a+-1*b*x), x) == :(-24*b^-1*exp(-1a+-1*b*x)+-1*b^-1*(a+b*x)^4*exp(-1a+-1*b*x)+-24*b^-1*(a+b*x)*exp(-1a+-1*b*x)+-12*b^-1*(a+b*x)^2*exp(-1a+-1*b*x)+-4*b^-1*(a+b*x)^3*exp(-1a+-1*b*x))
@test integrate((a+b*x)^4*(c+d*x)^-1*exp(-1a+-1*b*x), x) == :(-6*d^-1*exp(-1a+-1*b*x)+d^-4*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+d^-2*(-2*a*d+2*b*c)*exp(-1a+-1*b*x)+-1*d^-1*(a+b*x)^3*exp(-1a+-1*b*x)+-1*d^-3*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-6*d^-1*(a+b*x)*exp(-1a+-1*b*x)+-3*d^-1*(a+b*x)^2*exp(-1a+-1*b*x)+d^-5*(b*c+-1*a*d)^4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+d^-2*(a+b*x)^2*(b*c+-1*a*d)*exp(-1a+-1*b*x)+d^-2*(a+b*x)*(-2*a*d+2*b*c)*exp(-1a+-1*b*x)+-1*d^-3*(b*c+-1*a*d)^2*(a+b*x)*exp(-1a+-1*b*x))
@test integrate((a+b*x)^4*(c+d*x)^-2*exp(-1a+-1*b*x), x) == :(-2*b*d^-2*exp(-1a+-1*b*x)+-1*b^3*d^-4*(c+d*x)^2*exp(-1a+-1*b*x)+-1*d^-5*(c+d*x)^-1*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+-6*b*d^-4*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-2*b^2*d^-3*(c+d*x)*exp(-1a+-1*b*x)+4*b*d^-3*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-1*b*d^-6*(b*c+-1*a*d)^4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-4*b*d^-5*(b*c+-1*a*d)^3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+4*b^2*d^-4*(c+d*x)*(b*c+-1*a*d)*exp(-1a+-1*b*x))
@test integrate((a+b*x)^4*(c+d*x)^-3*exp(-1a+-1*b*x), x) == :(-1*b^2*d^-3*exp(-1a+-1*b*x)+b^2*d^-4*(-4*a*d+3*b*c)*exp(-1a+-1*b*x)+-1*x*b^3*d^-3*exp(-1a+-1*b*x)+-1//2*d^-5*(c+d*x)^-2*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+(1/2)*b*d^-6*(c+d*x)^-1*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+(1/2)*b^2*d^-7*(b*c+-1*a*d)^4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+4*b*d^-5*(c+d*x)^-1*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+4*b^2*d^-6*(b*c+-1*a*d)^3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+6*b^2*d^-5*(b*c+-1*a*d)^2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((a+b*x)^4*(c+d*x)^-4*exp(-1a+-1*b*x), x) == :(-1*b^3*d^-4*exp(-1a+-1*b*x)+-1//3*d^-5*(c+d*x)^-3*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+-6*b^2*d^-5*(c+d*x)^-1*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-6*b^3*d^-6*(b*c+-1*a*d)^2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-4*b^3*d^-5*(b*c+-1*a*d)*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-2*b^2*d^-6*(c+d*x)^-1*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+-2*b^3*d^-7*(b*c+-1*a*d)^3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+2*b*d^-5*(c+d*x)^-2*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+-1//6*b^2*d^-7*(c+d*x)^-1*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+-1//6*b^3*d^-8*(b*c+-1*a*d)^4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+1//6*b*d^-6*(c+d*x)^-2*(b*c+-1*a*d)^4*exp(-1a+-1*b*x))
@test integrate((a+b*x)^4*(c+d*x)^-5*exp(-1a+-1*b*x), x) == :(b^4*d^-5*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1//4*d^-5*(c+d*x)^-4*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+-3*b^2*d^-5*(c+d*x)^-2*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+3*b^3*d^-6*(c+d*x)^-1*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+3*b^4*d^-7*(b*c+-1*a*d)^2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+4*b^3*d^-5*(c+d*x)^-1*(b*c+-1*a*d)*exp(-1a+-1*b*x)+4*b^4*d^-6*(b*c+-1*a*d)*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-2//3*b^2*d^-6*(c+d*x)^-2*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+-1//24*b^2*d^-7*(c+d*x)^-2*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+1//12*b*d^-6*(c+d*x)^-3*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+1//24*b^3*d^-8*(c+d*x)^-1*(b*c+-1*a*d)^4*exp(-1a+-1*b*x)+1//24*b^4*d^-9*(b*c+-1*a*d)^4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+2//3*b^3*d^-7*(c+d*x)^-1*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+2//3*b^4*d^-8*(b*c+-1*a*d)^3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+4//3*b*d^-5*(c+d*x)^-3*(b*c+-1*a*d)^3*exp(-1a+-1*b*x))
@test integrate(F^(c*(a+b*x))*x^m*log(d*x)^n*(e+e*n+e*(1+m+b*c*x*log(F))*log(d*x)), x) == :(e*F^(c*(a+b*x))*x^(1+m)*log(d*x)^(1+n))
@test integrate(F^(c*(a+b*x))*x^2*log(d*x)^n*(e+e*n+e*(3+b*c*x*log(F))*log(d*x)), x) == :(e*F^(c*(a+b*x))*x^3*log(d*x)^(1+n))
@test integrate(x*F^(c*(a+b*x))*log(d*x)^n*(e+e*n+e*(2+b*c*x*log(F))*log(d*x)), x) == :(e*F^(c*(a+b*x))*x^2*log(d*x)^(1+n))
@test integrate(F^(c*(a+b*x))*log(d*x)^n*(e+e*n+e*(1+b*c*x*log(F))*log(d*x)), x) == :(e*x*F^(c*(a+b*x))*log(d*x)^(1+n))
@test integrate(F^(c*(a+b*x))*x^-1*log(d*x)^n*(e+e*n+b*c*e*x*log(F)*log(d*x)), x) == :(e*F^(c*(a+b*x))*log(d*x)^(1+n))
@test integrate(F^(c*(a+b*x))*x^-2*log(d*x)^n*(e+e*n+e*(-1+b*c*x*log(F))*log(d*x)), x) == :(e*F^(c*(a+b*x))*x^-1*log(d*x)^(1+n))
@test integrate(F^(c*(a+b*x))*x^-3*log(d*x)^n*(e+e*n+e*(-2+b*c*x*log(F))*log(d*x)), x) == :(e*F^(c*(a+b*x))*x^-2*log(d*x)^(1+n))
@test integrate(x^4*exp(a+b*x)^(1/2), x) == :(768*b^-5*exp(a+b*x)^(1/2)+-384*x*b^-4*exp(a+b*x)^(1/2)+-16*b^-2*x^3*exp(a+b*x)^(1/2)+2*b^-1*x^4*exp(a+b*x)^(1/2)+96*b^-3*x^2*exp(a+b*x)^(1/2))
@test integrate(x^3*exp(a+b*x)^(1/2), x) == :(-96*b^-4*exp(a+b*x)^(1/2)+-12*b^-2*x^2*exp(a+b*x)^(1/2)+2*b^-1*x^3*exp(a+b*x)^(1/2)+48*x*b^-3*exp(a+b*x)^(1/2))
@test integrate(x^2*exp(a+b*x)^(1/2), x) == :(16*b^-3*exp(a+b*x)^(1/2)+-8*x*b^-2*exp(a+b*x)^(1/2)+2*b^-1*x^2*exp(a+b*x)^(1/2))
@test integrate(x*exp(a+b*x)^(1/2), x) == :(-4*b^-2*exp(a+b*x)^(1/2)+2*x*b^-1*exp(a+b*x)^(1/2))
@test integrate(exp(a+b*x)^(1/2), x) == :(2*b^-1*exp(a+b*x)^(1/2))
@test integrate(x^-1*exp(a+b*x)^(1/2), x) == :(exp(a+b*x)^(1/2)*exp(-1//2*b*x)*expintegralei((1/2)*b*x))
@test integrate(x^-2*exp(a+b*x)^(1/2), x) == :(-1*x^-1*exp(a+b*x)^(1/2)+(1/2)*b*exp(a+b*x)^(1/2)*exp(-1//2*b*x)*expintegralei((1/2)*b*x))
@test integrate(x^-3*exp(a+b*x)^(1/2), x) == :(-1//2*x^-2*exp(a+b*x)^(1/2)+-1//4*b*x^-1*exp(a+b*x)^(1/2)+1//8*b^2*exp(a+b*x)^(1/2)*exp(-1//2*b*x)*expintegralei((1/2)*b*x))
@test integrate(x^-4*exp(a+b*x)^(1/2), x) == :(-1//3*x^-3*exp(a+b*x)^(1/2)+-1//12*b*x^-2*exp(a+b*x)^(1/2)+-1//24*b^2*x^-1*exp(a+b*x)^(1/2)+1//48*b^3*exp(a+b*x)^(1/2)*exp(-1//2*b*x)*expintegralei((1/2)*b*x))

max_stars_repo_path: test/2 Exponentials/2.1 u (F^(c (a+b x)))^n.jl


text:
using ModelingToolkitStandardLibrary.Magnetic, ModelingToolkit, OrdinaryDiffEq, Test

import ModelingToolkitStandardLibrary.Electrical
import ModelingToolkitStandardLibrary.Blocks
import ModelingToolkitStandardLibrary.Magnetic
using ModelingToolkit, OrdinaryDiffEq, Test
# using Plots

@parameters t

@testset "Inductor" begin
    mu_air = 1
    l_air = 0.0001
    mu_Fe = 1000
    l_Fe = 4*0.065
    a = b = 0.25

    @named source = Blocks.Sine(amplitude=230*sqrt(2), frequency=50, phase=pi/2)
    @named voltage = Electrical.Voltage()
    @named r = Electrical.Resistor(R=7.5)
    @named ground = Electrical.Ground()
    @named coil = Magnetic.FluxTubes.ElectroMagneticConverter(N=600)
    @named ground_m = Magnetic.FluxTubes.Ground()
    @named r_mAirPar = Magnetic.FluxTubes.ConstantReluctance(R_m=a * b * l_air * mu_air)
    @named r_mFe = Magnetic.FluxTubes.ConstantReluctance(R_m=a * b * l_Fe * mu_Fe)
    @named r_mLeak = Magnetic.FluxTubes.ConstantReluctance(R_m=1.2e6)
    connections = [
        connect(source.output, voltage.V)
        connect(voltage.p, r.p)
        connect(r.n, coil.p)
        connect(voltage.n, coil.n)
        connect(coil.port_p, r_mLeak.port_p)
        connect(r_mLeak.port_p, r_mAirPar.port_p)
        connect(r_mAirPar.port_n, r_mFe.port_p)
        connect(r_mFe.port_n, r_mLeak.port_n)
        connect(r_mFe.port_n, coil.port_n)
        connect(ground.g, voltage.n)
        connect(ground_m.port, r_mFe.port_n)
    ]
    @named model = ODESystem(connections, t, systems=[source, r, ground, coil, ground_m, r_mAirPar, r_mFe, r_mLeak, voltage])
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[], (0, 0.1))
    sol = solve(prob, Rodas4())

    # Plots.plot(sol; vars=[r.i])
    # Plots.plot(sol; vars=[r_mFe.V_m, r_mFe.Phi])

    @test sol.retcode == :Success
    @test sol[r_mFe.Phi] == sol[r_mAirPar.Phi]
    @test all(sol[coil.port_p.Phi] + sol[r_mLeak.Phi] + sol[r_mAirPar.Phi] .== 0)
end


text:
function getTestDataLinearInts(dataLength)
    firstN = Array{Int64,1}(undef, dataLength)
    for ii in 1:dataLength
        firstN[ii] = ii
    end
    return firstN
end
function getTestDataRandomInts(dataLength)
    firstNrand = Array{Int64,1}(undef, dataLength+1)
    randDiff = Array{Int64,1}(undef, dataLength)
    firstNrand[1] = rand(1:dataLength)
    for ii in 1:dataLength
        firstNrand[ii+1] = rand(1:dataLength)
        randDiff[ii] = firstNrand[ii+1] - firstNrand[ii]
    end
    return randDiff
end
function getTestDataLinearFloats(dataLength)
    firstN = Array{Float32,1}(undef, dataLength)
    startN = 5.0
    lastN = 10.0
    dataRange = lastN - startN;
    dataStep = dataRange/(dataLength-1)
    firstN[1] = startN
    for ii in 2:dataLength
        firstN[ii] = firstN[ii-1] + dataStep
    end
    return firstN
end


text:
include("../tatg.jl")
using CompEcon: qnwnorm
using JLD
using KernelDensity
using PyCall
using PyPlot
using StatsBase
@pyimport seaborn as sb

# ------------------------------------------------------------------- #
# File Helper Functions
# ------------------------------------------------------------------- #
"Reads lj, lu, lλ☆p, m_ls from a specific file"
function read_solutions(filename::AbstractString)
    file = jldopen(filename, "r")
    lj, lu, lλ☆p, m_ls =
        map(var->read(file, var), ["lj", "lu", "llambdastar", "m_ls"]);
    close(file)

    return lj, lu, lλ☆p, m_ls
end

"""
Takes model parameters and a filename and returns a model with that has
stochastic volatility with the correct parmeters and all relevant
solution output
"""
function build_modelSV_interpolants(α::Float64, ρ::Float64, σ::Float64;
                                    file="solutions/curr_soln_sv_9.jld")
    # Get solutions
    lj, lu, lλ☆p, m_ls = read_solutions(file)

    # Get inputs to model creation
    nλ☆ = length(m_ls[1])
    nlzh = length(m_ls[2])
    nv = length(m_ls[3])
    lλ☆bds = m_ls[1][end]
    neps = Int(size(lλ☆p, 1)^(1/3))

    # Create model
    m = Exchange.ModelSV(α1=α, α2=α, lλ☆bds=lλ☆bds,
                         nλ☆=nλ☆, nlzh=nlzh, nv=nv,
                         nϵ1=neps, nϵ2=neps, nϵ3=neps)

    # Get interpolants
    vf_itp = Exchange.update_vf_itps(m, lj, lu);
    lλ☆p_itp = Exchange.get_lλ☆p_itp(m, lλ☆p);

    return m, vf_itp, lλ☆p_itp, lj, lu, lλ☆p
end

"""
Takes model parameters and a filename and returns a model with that has
constant volatility with the correct parmeters and all relevant
solution output
"""
function build_modelCV_interpolants(α::Float64, ρ::Float64, σ::Float64;
                                    filename="solutions/curr_soln_cv_9.jld")
    # Get solutions
    lj, lu, lλ☆p, m_ls = read_solutions(filename)

    # Get inputs to model creation
    nλ☆ = length(m_ls[1])
    nlzh = length(m_ls[2])
    lλ☆bds = m_ls[1][end]
    neps = Int(size(lλ☆p, 1)^(1/2))

    # Create model
    m = Exchange.ModelCV(α1=α, α2=α, lλ☆bds=lλ☆bds,
                         nλ☆=nλ☆, nlzh=nlzh,
                         nϵ1=neps, nϵ2=neps)

    # Get interpolants
    vf_itp = Exchange.update_vf_itps(m, lj, lu);
    lλ☆p_itp = Exchange.get_lλ☆p_itp(m, lλ☆p);

    return m, vf_itp, lλ☆p_itp, lj, lu, lλ☆p
end

# ------------------------------------------------------------------- #
# General Helper Functions
# ------------------------------------------------------------------- #
"""
Simulates for a long time with no shocks in order to go towards the
steady state
"""
function find_steadystate_lλ(m::Exchange.AbstractModel,
                             lλ☆p_itp::Exchange.AbstractInterpolation;
                             capT=2_000_000)

    return Exchange.simulate(m, lλ☆p_itp, zeros(3, capT))[7][end]
end

"""
Takes the ratio of consumptions, an exchange rate, and an alpha and returns
slope and correlations that are implied
"""
function get_slope_corr(c2_c1_ratio, fxr, alpha, npers=5000)
    x = log(c2_c1_ratio[1:npers])
    y = log(fxr[1:npers])
    slope = x \ y
    corrr = round(cor(log(c2_c1_ratio[1:npers]), log(fxr[1:npers])), 2)

    return slope, corrr
end

# ------------------------------------------------------------------- #
# Pareto Frontiers
# ------------------------------------------------------------------- #
"""
Plots the Pareto frontiers for both consumption and value functions
"""
function plot_c_vf_frontiers!(lj, lu, c1, c2, ax)

    # Set colors and line style
    colors = ("#852C64", "#AA6039")
    ls = ("-", "-.")

    # Get how many of the exogenous states are and find mean
    i = map(x->div(x, 2)+1, size(lj)[2:end])

    ax[:plot](exp(lu[:, i...]), exp(lj[:, i...]), color=colors[1], ls=ls[1])
    ax[:plot](c2[:, i...], c1[:, i...], color=colors[2], ls=ls[2])

    ax[:set_ylabel]("Agent 1")
    ax[:set_xlabel]("Agent 2")
    return ax
end

"""
Plots the frontier for value functions only
"""
function plot_vf_frontiers!(lj, lu, c, ax)

    ax[:plot](exp(lu), exp(lj), color=c)
    ax[:set_xlabel](L"log U")
    ax[:set_ylabel](L"log J")

    return ax
end

# ------------------------------------------------------------------- #
# Consumption v Exchange Rate Functions
# ------------------------------------------------------------------- #
"""
Takes a vectors of `c2_c1_ratio` and a vector of `fxr` and plots them
on a scatter plot.
"""
function plot_c2c1ratio_fxr_single!(alpha, c2_c1_ratio, fxr, npers, ax)
    # Get slope and correlation
    s, c = get_slope_corr(c2_c1_ratio, fxr, alpha, npers);

    palpha = npers < 10000 ? .75 : npers < 25000 ? .5 : .25

    ax[:scatter](log(c2_c1_ratio[1:npers]), log(fxr[1:npers]);
                 s=0.75, alpha=1.00, label="Correlation is $c", plotkwargs[alpha]...)

    return ax
end

# ------------------------------------------------------------------- #
# Pareto Weight Plots
# ------------------------------------------------------------------- #
"""
Plots the history of Pareto weights. Has option of plotting a random
walk of the same standard deviation next to history.
"""
function plot_pareto_history!(alpha, lλ☆, ax; rw=false, T::Int=1_000_000,
                              stsz::Int=10000, c="red")

    ax[:plot](collect(1:stsz:T), lλ☆[1:stsz:T]; color=c)

    if rw
        srand(4251989)
        lλ☆_std = std(lλ☆[1:end-1] - lλ☆[2:end])
        ranwalk = cumsum(randn(T) .* lλ☆_std)
        ax[:plot](collect(1:stsz:T), ranwalk[1:stsz:T]; color="k")
    end

    ax[:set_ylim](-3.75, 3.75)

    return ax
end

# ------------------------------------------------------------------- #
# Exchange Rate
# ------------------------------------------------------------------- #
"""
Plots the autocorrelation function for a specified number of periods
"""
function acorr_plot_one_ax!(data, α, ax, k_max=100, nshow=30; c="red")
    ks = collect(1:k_max)
    acf = autocor(data, ks)
    lab = string(L"$\rho^k$ ($\rho=", round(acf[1], 4), L")$")
    ax[:plot](ks, acf; label=lab, color=c)
    ax
end

# ------------------------------------------------------------------- #
# Impulse Response Function
# ------------------------------------------------------------------- #
"""
Plots linear impulse response functions. In particular, given a history
of shocks, ϵ, it plots the response of (zhat, v, c1, c2, lλ☆, e) in %
deviations
"""
function irf_plotter(m::Exchange.AbstractModel,
                     lλ☆p_itp::Exchange.Interpolations.AbstractInterpolation,
                     ϵ::Array{Float64, 2};
                     capT=50, lλ0=0.0, c="red")

    # Get irf data
    simuls = Exchange.simulate(m, lλ☆p_itp, ϵ; lλ0=lλ0)
    a1, b1, c1, a2, b2, c2, lλ☆, exogvals, c2_c1_ratio, fxr = simuls

    # Create a t array
    irf_per = collect(1:capT)
    fig_irf, ax_irf = subplots(3, 2)
    map(x->x[:locator_params](nbins=4), ax_irf)

    # Put zhat in log deviations and plot
    ld_zhat = 100*(exogvals[1, 1:end-1]' - exogvals[1, 1])
    ax_irf[1, 1][:plot](irf_per, ld_zhat; color="k")
    ax_irf[1, 1][:set_title]("Productivity")

    # Put v in log deviations and plot
    ld_v = 100*(log(exogvals[2, 1:end-1]') - log(exogvals[2, 1]))
    ax_irf[1, 2][:plot](irf_per, ld_v, color="k")
    ax_irf[1, 2][:set_title]("Volatility in Country 1")

    # Put c1 in log deviations and plot
    ld_c1 = 100*(log(c1) - log(c1[1]))
    ax_irf[2, 1][:plot](irf_per, ld_c1; color=c)
    ax_irf[2, 1][:set_title]("Consumption: Country 1")

    # Put c2 in log deviations and plot
    ld_c2 = 100*(log(c2) - log(c2[1]))
    ax_irf[2, 2][:plot](irf_per, ld_c2; color=c)
    ax_irf[2, 2][:set_title]("Consumption: Country 2")

    # Put P.W. ratio in log deviations and plot
    ld_λ☆ = 100*(lλ☆ - lλ☆[1])
    ax_irf[3, 1][:plot](irf_per, ld_λ☆; color=c)
    ax_irf[3, 1][:set_title]("Pareto Weight Ratio")
    ax_irf[3, 1][:set_xlabel]("Time (Quarters)")

    # Put exchange rate in log deviations and plot
    ld_fxr = 100*(log(fxr) - log(fxr[1]))
    ax_irf[3, 2][:plot](irf_per, ld_fxr; color=c)
    ax_irf[3, 2][:set_title]("Exchange Rate")
    ax_irf[3, 2][:set_xlabel]("Time (Quarters)")

    fig_irf[:tight_layout]()
    # fig_irf[:savefig]("images/irf_wrt_v.png", dpi=400)
    fig_irf[:show]()

    return fig_irf, ax_irf
end

# ------------------------------------------------------------------- #
# Policies
# ------------------------------------------------------------------- #
"""
Plots the policy function over the possible values for vt. All other
variables are held at their ss.
"""
function plot_policy_wrt_v!(m, itp, ax=nothing; c="red")
    lλ☆t = 0.0
    lzh = 0.0
    v = linspace(extrema(m.grid[:, 3])..., 50)

    lλ☆p = Array(Float64, 50)

    for i=1:50
        lλ☆p[i] = itp[lλ☆t, lzh, v[i], zeros(3)...]
    end

    if ax === nothing
        fig, ax = subplots()
    end

    ax[:set_ylabel](L"$\log \lambda_{t+1}^*$")
    ax[:set_xlabel](L"$v_{t}$")
    ax[:plot](v, lλ☆p; color=c)
    ax
end

"""
Plots the expected change in the Pareto weight across various values of
lλ☆. All other variables are held at their steady states.
"""
function plot_policy_at_ss!(m, itp, ss, ax=nothing; c="red")

    lλ☆t = linspace(ss-3.5, ss+3.5, 50)
    lzh = 0.0
    v = m.exog.vbar

    ϵ, Π = Exchange.qnwnorm([7, 7, 7], zeros(3), eye(3))
    nshocks = length(Π)
    lλ☆p = Array(Float64, nshocks, 50)

    for i=1:50
        for j=1:nshocks
            lλ☆p[j, i] = itp[lλ☆t[i], lzh, v, ϵ[j, :]...]
        end
    end

    Elλ☆p = vec(Π' * lλ☆p)
    ΔElλ☆p = Elλ☆p - collect(lλ☆t)

    if ax === nothing
        fig, ax = subplots()
    end

    ax[:set_ylabel](L"$E[\log \lambda_{t+1}^*] - \log \lambda_{t}^*$")
    ax[:set_xlabel](L"$\log \lambda_{t}^*$")
    ax[:plot](lλ☆t, ΔElλ☆p; color=c)
    ax
end


text:
using Test
using QuantumLattices.Mathematics.Combinatorics

@testset "Combinations" begin
    @test Combinations{0}(1:3)|>collect == [()]
    @test Combinations{1}(1:3)|>collect == [(1,), (2,), (3,)]
    @test Combinations{2}(1:3)|>collect == [(1, 2), (1, 3), (2, 3)]
    @test Combinations{3}(1:3)|>collect == [(1, 2, 3)]
end

@testset "DulCombinations" begin
    @test DulCombinations{0}(1:3)|>collect == [()]
    @test DulCombinations{1}(1:3)|>collect == [(1,), (2,), (3,)]
    @test DulCombinations{2}(1:3)|>collect == [(1, 1), (1, 2), (1, 3), (2, 2), (2, 3), (3, 3)]
    @test DulCombinations{3}(1:3)|>collect == [(1, 1, 1), (1, 1, 2), (1, 1, 3), (1, 2, 2), (1, 2, 3), (1, 3, 3), (2, 2, 2), (2, 2, 3), (2, 3, 3), (3, 3, 3)]
end

@testset "Permutations" begin
    @test Permutations{0}(1:3)|>collect == [()]
    @test Permutations{1}(1:3)|>collect == [(1,), (2,), (3,)]
    @test Permutations{2}(1:3)|>collect == [(1, 2), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2)]
    @test Permutations{3}(1:3)|>collect == [(1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1)]
end

@testset "DulPermutations" begin
    @test DulPermutations{0}(1:3)|>collect == [()]
    @test DulPermutations{1}(1:3)|>collect == [(1,), (2,), (3,)]
    @test DulPermutations{2}(1:3)|>collect == [(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3)]
    @test DulPermutations{3}(1:3)|>collect == [(i, j, k) for i = 1:3 for j = 1:3 for k = 1:3]
end


text:
using Plots
L=100
p=0.8
function LimitControl(arr, i, j)
    color_change_list=[]#this list includes the neighbors' attributes
    if i+1!=L+1
        push!(color_change_list,arr[i+1,j])
    end
    if i-1!=0
        push!(color_change_list,arr[i-1,j])
    end
    if j+1!=L+1
        push!(color_change_list,arr[i,j+1])
    end
    if j-1!=0
        push!(color_change_list, arr[i,j-1])
    end
    return color_change_list
end
function whole_color_change(arr,color_min, row, col, color_change_list)
    arr[row,col]=color_min
    for i in 1:L
        for j in 1:L
            for color in color_change_list
                if arr[i,j]==color
                    arr[i,j]=color_min
                end
            end
        end
    end
    return arr
end
function neighbor_checking(arr, row, col,color_change_list)
    color_min=minimum(color_change_list)
    count=0
    the_neighbor_color=0
    for i in color_change_list
        if i==0
            return arr
            break
        end
        if i!=0
            count+=1
            the_neighbor_color=i
        end
    end
    if count==1 #if true, the_neighbor_color is the color of that only neighbor.
        arr[row, col]=the_neighbor_color
        return  arr
    end
    if count>1
        arr=whole_color_change(arr,color_min, row, col, color_change_list)
    end
        return arr
end
function Percolation(L,p)
    tot=[]
    int_max=100000
    color_=2
    arr=zeros((L,L))
    for i in 1:L
        arr[i,1]=1
        arr[i,L]=int_max
    end
    for j in 2:L-1
        for i in 1:L
            color_change_list=LimitControl(arr, i , j)
            if p>rand()
                arr[i,j]=color_
                color_+=1
            end
            arr=neighbor_checking(arr, i, j,color_change_list) #color_change_list includes the neighbors' color, so we can use the list as the neighbors' attributes.
        end
        push!(tot, copy(arr))
    end
    return tot
end
tot=Percolation(L,p)

@gif for i in 1:1:length(tot)
    heatmap(tot[i][:,1:end-1], c = cgrad(:copper, 10), legend = false, border = :none, title = "Percolation for P = 0.5")
end


text:
function generate_component!(block, module_name, prefix, meta)
    args = isempty(meta["args"]) ? Symbol[] : Symbol.(meta["args"])
    wild_args = isempty(meta["wild_args"]) ? Symbol[] : Symbol.(meta["wild_args"])
    fname = string(prefix, "_", lowercase(meta["name"]))
    fsymbol = Symbol(fname)

    append!(block.args,
            (quote
                export $fsymbol
                function $(fsymbol)(;kwargs...)
                    available_props = $args
                    wild_props = $wild_args
                    return Component($fname, $(meta["name"]), $module_name, available_props, wild_props; kwargs...)
                end
            end).args
    )
    signatures = String[string(repeat(" ", 4), fname, "(;kwargs...)")]
    if in(:children, args)
        append!(block.args,
            (quote
                $(fsymbol)(children::Any; kwargs...) = $(fsymbol)(;kwargs..., children = children)
                $(fsymbol)(children_maker::Function; kwargs...) = $(fsymbol)(children_maker();kwargs...)
            end).args
        )
        push!(signatures,
            string(repeat(" ", 4), fname, "(children::Any, kwargs...)")
        )
        push!(signatures,
            string(repeat(" ", 4), fname, "(children_maker::Function, kwargs...)")
        )
    end
    docstr = string(
        join(signatures, "\n"),
        "\n\n",
        meta["docstr"]
    )
    push!(block.args, :(@doc $docstr $fsymbol))
end

function generate_components_package(meta)
    result = Expr(:block)
    name  = meta["name"]
    prefix  = meta["prefix"]
    for cmeta in meta["components"]
        generate_component!(result, name, prefix, cmeta)
    end
    return result
end

function generate_embeded_components()
    dash_meta = load_meta("dash")
    packages = dash_meta["embedded_components"]
    result = Expr(:block)
    for p in packages
        append!(result.args,
            generate_components_package(load_meta(p)).args
        )
    end
    return result
end

macro place_embedded_components()
    return esc(
        generate_embeded_components()
    )
end


text:
import TimeSeries: TimeArray

function varEstimate(data::TimeArray{Float64}, lags::Int, typ::AbstractString)
	t = varEstimate(data.values, lags, typ)
	t.seriesNames = data.colnames
	return t
end


text:
#!/usr/bin/env julia
#
# Benchmarking script for singular value problems

try
    Pkg.installed("Benchmarks")
catch
    Pkg.clone("https://github.com/johnmyleswhite/Benchmarks.jl")
end

try
    Pkg.installed("PROPACK")
catch
    Pkg.clone("https://github.com/andreasnoack/PROPACK.jl")
end

using Benchmarks
using IterativeSolvers
using JLD
using MAT
using PROPACK

BASEDIR = "florida"

if !isdir(BASEDIR)
    info("Setting up $BASEDIR collection")
    include("setup-florida.jl")
end

"""
Run benchmark function `runbenchmark` against all .mat files in all
subdirectories of BASEDIR
"""
function benchmarkall(BASEDIR)
    for group in readdir(BASEDIR)
        isdir(joinpath(BASEDIR, group)) || continue
        for matrix in readdir(joinpath(BASEDIR, group))
            filename = joinpath(BASEDIR, group, matrix)

            endswith(matrix, ".mat") || continue

        #If we already saved benchmarks, don't rerun
        benchmarkfilename = joinpath(BASEDIR, group, matrix[1:end-3]*"jld")
        isfile(benchmarkfilename) && continue

        #To debug this script, it's useful to run it on small matrices only
        #filesize(filename) < 100000 || continue

        runbenchmark(filename, benchmarkfilename)

        #debug: just try one
        #@goto stop
        end
    end
    @label stop
end

"""
Run several different computations for singular values
"""
function runbenchmark(filename, benchmarkfilename)
    mf = matread(filename)
    if !haskey(mf, "Problem")
        warn("Skipping unknown file $filename: No 'Problem' struct")
        return
    end

    pr = mf["Problem"]
    if !haskey(pr, "A")
        warn("Skipping unknown file $filename: 'Problem' struct has no matrix 'A'")
        return
    end
    A = pr["A"]
    #eltype(A) <: Real || warn("Skipping matrix with unsupported element type $(eltype(A))")
    info(filename*", matrix of type $(typeof(A)) and size $(size(A))")
    m, n = size(A)
    #Choose the same normalized unit vector to start with
    q = randn(n)
    eltype(A) <: Complex && (q += im*randn(n))
    rmul!(q, inv(norm(q)))
    qm = randn(m)
    eltype(A) <: Complex && (q += im*randn(m))
    rmul!(qm, inv(norm(qm)))

    #Number of singular values to request
    nv = min(m, n, 10)

    #Maximum number of iterations
    maxiter = max(m, n)

    #Tolerance, however the algorithm chooses to interpret this
    #Set convergence criterion to sqrt eps
    tol = √eps(real(one(eltype(A))))

    @info("svds (eigs on [0 A; A' 0])")
    #See #12890 :(
    #b_svds = @benchmark svds(A, v0=q, nsv=nv, tol=tol, maxiter=maxiter)
    b_svds = @benchmark svds(A, nsv=nv, tol=tol, maxiter=maxiter)

    @info("eigs on A'A or AA'")
    makeata(A) = (size(A, 1)≥ size(A,2) ? A*A' : A'A)
    B = makeata(A)
    b_ata = @benchmark makeata(A)
    b_eigs = if m≥n
        B = A'A #Size n x n
        b_eigs = @benchmark eigs(B, v0=q, nev=nv, tol=tol, maxiter=maxiter)
    else
        B = A*A' #Size m x m
        #eigs cannot use the same starting vector as everyone else, needs R^m not R^n
        b_eigs = @benchmark eigs(B, v0=qm, nev=nv, tol=tol, maxiter=maxiter)
    end

    @info("tsvd (PROPACK)")
    b_tsvd = try
        if m≥n #PROPACK implicitly assumes the input is tall
            @benchmark tsvdvals(A, initvec=qm, kmax=maxiter, k=nv, tolin=tol)
        else
            @benchmark tsvdvals(A', initvec=q, kmax=maxiter, k=nv, tolin=tol)
        end
    catch exc
        println("Exception: $exc")
    end

    @info("GKL with thick restart using Ritz values")
    b_tr = try
        @benchmark svdl(A, nv, v0=q, tol=tol, reltol=tol)
    catch exc
        println("Exception: $exc")
    end

    @info("GKL with thick restart using harmonic Ritz values")
        b_trh = try
            @benchmark svdl(A, nv, v0=q, tol=tol, reltol=tol, method=:harmonic)
        catch exc
            println("Exception: $exc")
        end

    data = Dict(
        "svds" => b_svds,
        "ata"  => b_ata,
        "eigs" => b_eigs,
        "tsvd" => b_tsvd,
        "tr"   => b_tr,
        "trh"  => b_trh
    )
    println("Timings:")
    for (k, v) in data
        println(k)
        display(v)
        println()
    end
    JLD.save(benchmarkfilename, data, compress=true)
end

benchmarkall(BASEDIR)


text:
"""
Subroutine dicfs

Given a symmetric matrix A in compreessed column storage, this
subroutine computes an incomplete Cholesky factor of A + alpha*D,
where alpha is a shift and D is the diagonal matrix with entries
set to the l2 norms of the columns of A.

MINPACK-2 Project. October 1998.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function dicfs(n, nnz, A, L,
               p, alpha, iwa, wa1, wa2)
    nbmax = 3
    alpham = 1.0e-3
    nbfactor = 512

    zero = 0.0
    one = 1.0
    two = 2.0

    # Compute the l2 norms of the columns of A.
    nrm2!(wa1, A, n)

    # Compute the scaling matrix D.
    @inbounds for i=1:n
        if wa1[i] > zero
            wa2[i] = one/sqrt(wa1[i])
        else
            wa2[i] = one
        end
    end

    # Determine a lower bound for the step.

    if alpha <= zero
        alphas = alpham
    else
        alphas = alpha
    end

    # Compute the initial shift.

    alpha = zero
    @inbounds for i=1:n
        if getdiagvalue(A,i) == zero
            alpha = alphas
        else
            alpha = max(alpha,-getdiagvalue(A,i)*(wa2[i]^2))
        end
    end

    if alpha > 0
        alpha = max(alpha,alphas)
    end

    # Search for an acceptable shift. During the search we decrease
    # the lower bound alphas until we determine a lower bound that
    # is not acceptaable. We then increase the shift.
    # The lower bound is decreased by nbfactor at most nbmax times.

    nb = 1
    while true

        if isa(A, TronSparseMatrixCSC)
            # Copy the sparsity structure of A into L.
            @inbounds for i=1:n+1
                L.colptr[i] = A.colptr[i]
            end
            @inbounds for i=1:nnz
                L.rowval[i] = A.rowval[i]
            end

            # Scale A and store in the lower triangular matrix L.
            @inbounds for j=1:n
                L.diag_vals[j] = A.diag_vals[j]*(wa2[j]^2) + alpha
            end
            @inbounds for j=1:n
                for i=A.colptr[j]:A.colptr[j+1]-1
                    L.tril_vals[i] = A.tril_vals[i]*wa2[j]*wa2[A.rowval[i]]
                end
            end
        else
            # TronDenseMatrix case.
            @inbounds for j=1:n
                L.vals[j,j] = A.vals[j,j]*(wa2[j]^2) + alpha
            end
            @inbounds for j=1:n,k=j+1:n
                L.vals[k,j] = A.vals[k,j]*wa2[j]*wa2[k]
            end
        end

        # Attempt an incomplete factorization.
        info = dicf(n,nnz,L, p,
                   view(iwa,1:n),view(iwa,n+1:2*n),
                   view(iwa,2*n+1:3*n),wa1)

        # If the factorization exists, then test for termination.
        # Otherwise increment the shift.
        if info >= 0
            # If the shift is at the lower bound, reduce the shift.
            # Otherwise undo the scaling of L and exit.
            if alpha == alphas && nb < nbmax
                alphas /= nbfactor
                alpha = alphas
                nb = nb + 1
            else
                if isa(L, TronSparseMatrixCSC)
                    @inbounds for i=1:n
                        L.diag_vals[i] /= wa2[i]
                    end
                    @inbounds for j=1:L.colptr[n+1]-1
                        L.tril_vals[j] = L.tril_vals[j]/wa2[L.rowval[j]]
                    end
                else
                    @inbounds for j=1:n,k=j:n
                        L.vals[k,j] /= wa2[k]
                    end
                end
                return
            end
        else
            alpha = max(two*alpha,alphas)
        end
    end

    return
end


text:
module Ballistics

# Write your package code here.
#1. Exports/imports
export Fragmentation,Types, ExternalBallistics, Vulnerability, scenario, materials, Utils, inputs
#include("inputs.jl")
#using .inputs
include("Types.jl")
using .Types
include("materials.jl")
using .materials
include("scenario.jl")
using .scenario
include("Utils.jl")
using .Utils
#include("aerodynamic_coeff30.jl")
#using .aero30
#include("aerodynamic_coeff105.jl")
#using .aero105
include("InternalBallistics.jl")
using .InternalBallistics
include("ExternalBallistics.jl")
using .ExternalBallistics
include("fragmentation.jl")
using .Fragmentation
include("Vulnerability.jl")
using .Vulnerability
include("post.jl")
using .post

#2. Interface documentation
# in ballistics computation we have a projectile and a target
# a projectile accelerate with the
# InternalBallistcs!(p) function and flies with the
# ExternalBallistics!(p) function  and perfom penetration with
# Vulnerability!(p) function

#3. Generic definition for the interface
#Hard contracts
"""
InternalBallistics to get muzzle velocity. Optional
"""
function InternalBallistics!(args...)
nothing
end
"""
ExternalBallistics to get impact conditions. Optional
"""
function ExternalBallistics!(args...)
nothing
end
"""
Vulnerability to get damage assessment.
"""
function Vulnerability! end

#4. Game logic

function ballisticsComp!(p::Projectile1D, t::Target1D)
    InternalBallistics!(p)
    ExteralBallistics!(p)
    Vulnerability!(p,t)
    nothing
end
end


text:
using Test
@testset "QFTBlock" begin
    include("QFTBlock.jl")
end

@testset "RotBasis" begin
    include("RotBasis.jl")
end

@testset "Sequence" begin
    include("Sequence.jl")
end

@testset "Diff" begin
    include("Diff.jl")
end

@testset "Bag" begin
    include("Bag.jl")
end

@testset "TrivilGate" begin
    include("TrivilGate.jl")
end

@testset "ConditionBlock" begin
    include("ConditionBlock.jl")
end

@testset "Mod" begin
    include("Mod.jl")
end

@testset "pauli_strings" begin
    include("pauli_strings.jl")
end

@testset "reflect_gate" begin
    include("reflect_gate.jl")
end

@testset "math_gate" begin
    include("math_gate.jl")
end

@testset "shortcuts" begin
    include("shortcuts.jl")
end


text:
# module Poly

import Base: ==, !=, isless, isconst, isreal, isnull, isequal
import Base: +, -, *, /, ^, conj, conj!, abs2, norm, real, imag
import Base: show, print, convert, copy, hash, merge
import Base: start, next, done, length, setindex!, getindex, haskey, keys, values

abstract type  AbstractPolynomial end

"""
    Degree(explvar::Int, conjvar::Int)

Define a mathematical degree, which is used to define a variable exponent to
`explvar` and its conjugate exponent to `conjvar`.

### Attributes
- `explvar`: Int
- `conjvar`: Int
"""
mutable struct Degree
    explvar::Int
    conjvar::Int
end

"""
    Variable(varname::String, Complex)

Define a mathematical variable `varname` of a certain mathematical kind,
`Complex` here.

### Attributes
- `name`: String
- `kind`: a type among `Complex`, `Real` and `Bool`
"""
mutable struct Variable <: AbstractPolynomial
    name::String
    kind::Type

    function Variable(name, kind)
        if kind ∉ Set([Complex, Real, Bool]) || typeof(name) ∉ Set([String, SubString{String}])
            error("Variable() : attempting to define a variable $name of type $kind, supported types are {Complex, Real, Bool}")
        end
        return new(String(name), kind)
    end
end

"""
    Exponent(expo::Dict{Variable, Degree})

Define a mathematical exponent, that is a product of `Variable` and conjugated
`Variable`.

### Attributes
- `expo` : a dictionary associating `Variable` objects to a `Degree` objects,
their exponent and the exponent of their conjugates.
- `degree`: a `Degree` object indicating the global degree of the exponent :
the sum of the explicit variables degrees and of the conjugated variables
degrees.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > Exponent(Dict(a=>Degree(1,0), b=>Degree(1,1))) == a*b*conj(b)
```
"""
struct Exponent <: AbstractPolynomial
    expo::Dict{Variable, Degree}
    degree::Degree

    function Exponent(expo::Dict{Variable, Degree})
        degexpl, degconj = 0,0
        for (var, degree) in expo
            ((degree.explvar < 0) || (degree.conjvar < 0)) && error("Exponent(): Expected non negative exponent for variable $var (got $degree)")
            (isreal(var) && degree.conjvar != 0) && error("Exponent(): Expected nul conj exponent for real variable $var (got $degree)")
            (isbool(var) && degree.explvar ∉ Set([0,1])) && error("Exponent(): Expected boolean exponent for bool $var (got $degree)")
            if degree != Degree(0,0)
                degexpl += degree.explvar
                degconj += degree.conjvar
            else
                delete!(expo, var)
            end
        end
        return new(expo, Degree(degexpl, degconj))
    end
end

"""
    Polynomial(poly::Dict{Exponent, Number}, degree::Degree)

Define a mathematical polynomial, that is a linear combinaison of product of
`Variable` and conjugated `Variable`.

### Attributes
- `poly` : a dictionary associating an `Exponent` to a a number.
- `degree`: a `Degree` object indicating the global degree of the polynomial :
the maximum of the explicit exponent degrees and of the conjugated exponent
degrees.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > p = a*b*conj(b) + (2+3im) * b^3
julia > p.poly = Dict(a*b*conj(b)=>1, b^3=>2+3im)
julia > p.degree = (3,1)
```
"""
struct Polynomial <: AbstractPolynomial
    poly::Dict{Exponent, Number}
    degree::Degree

    function Polynomial(poly::Dict{Exponent, Number})
        degexpl, degconj = 0, 0
        for (expo, λ) in poly
            if λ!=0
                degexpl = max(degexpl, expo.degree.explvar)
                degconj = max(degconj, expo.degree.conjvar)
            else
                delete!(poly, expo)
            end
        end
        return new(poly, Degree(degexpl, degconj))
    end
end

"""
    Point(coords::Dict{Variable, Number})

Define a mathematical point, that is a pairing of variables and numbers.

### Attributes
- `coords` : a dictionary associating a `Variable` to a a number.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > pt = Point(Dict(a=>π, b=>e+7im))
```
"""
struct Point
    coords::Dict{Variable, Number}

    function Point(dict::Dict)
        dict_pt = Dict{Variable, Number}()
        for (var, val) in dict
            if !isa(var, Variable) || !isa(val, Number)
                error("Point(): Expected pair of (Variable, Number), got ($var, $val) of type ($typeof(var), $typeof(val)) instead.")
            end
            if isbool(var)
                booled = 0
                if val != 0
                    booled = Int((val/abs(val) + 1) / 2)
                end
                (booled ≠ val) && warn("Point(): $var is $(var.kind), provided value is $val, $booled stored.")
                add_to_dict!(dict_pt, var, booled)
            elseif isreal(var)
                realed = real(val)
                (realed ≠ val) && warn("Point(): $var is $(var.kind), provided value is $val, $realed stored.")
                add_to_dict!(dict_pt, var, realed)
            else
                add_to_dict!(dict_pt, var, val)
            end
        end
        return new(dict_pt)
    end
end

"""
    Constraint(p::AbstractPolynomial, lb::Number, ub::Number)

Define a mathematical complex constraint, from a polynomial body and two complex
bounds. Equivalent to the two real inequality produced from real and imag part
of `p`, `lb` and `ub`.

### Attributes
- `p::Polynomial` : the polynomial body of the constraint.
- `lb::Number` : complex lower bound.
- `ub::Number` : complex upper bound.
- `precond::Symbol` : a symbol describing the preconditioning method to be
applied for this constraint. Default value is `:none`.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > cstr = abs2(a^2) + abs2(3*b) << 25
julia > cstr.precond = :sqrt
```
"""
mutable struct Constraint
  p::Polynomial
  lb::Number
  ub::Number
  precond::Symbol

  Constraint(p::AbstractPolynomial, lb::Number, ub::Number) = new(p, lb, ub, :none)
end

"""
    Problem()

Define an empty mathematical polynomial optimization problem, characterized by
a polynomial `objective`, a dictionnary of `constraints` and a dictionnary of
used variables `variables`.

### Attributes
- `objective::Polynomial` : polynomial criterion.
- `constraints::Dict{String, Constraint}` : dictionnary of constraint name to
constraint.
- `variables::Dict{String, Type}` : dictionnary of variable name to type.

### Exemple
```julia
julia > x, y, z = ...
julia > pb = Problem()
julia > set_objective(pb, abs2(x+y+z))
julia > add_constraint!(pb, abs2(x) << 1)
julia > add_constraint!(pb, y << 0.5+1im)
julia > add_constraint!(pb, 0-1im << z << 1+0im)
```
"""
mutable struct Problem
  objective::Polynomial
  constraints::Dict{String, Constraint}
  variables::Dict{String, Type}
end

include("Modeler_accessors.jl")
include("Modeler_constructors.jl")
include("Modeler_cplx2real.jl")
include("Poly_Cplx2Real.jl")
include("Poly_operators.jl")
include("Poly_order.jl")
include("Poly_point_algebra.jl")
include("Poly_poly_algebra.jl")
include("Poly_print.jl")
include("utils_ampl.jl")
include("utils_dat_compare.jl")
include("utils_dat_export.jl")
include("utils_dat_import.jl")
include("utils_Poly.jl")

# export Variable, Point, Exponent, Monomial, Polynomial
# export isconst, isone, is_homogeneous
# export evaluate, abs2, conj, add!, add
# export cplx2real, real, imag
# export norm

#end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using InvertibleNetworks, NNlib, LinearAlgebra, Test

# Data
nx = 28
ny = 28
n_in = 8
batchsize = 1
n_hidden = 8
k = 3   # kernel size
s = 1   # stride
p = 1   # padding

pars = (nx, ny, n_in, batchsize, k, s,p)
###################################################################################################
# Test layer invertibility

function test_inv(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    print("\nInvertibility test hyperbolic layer with action=$(action)\n")
    # Data
    X0 = randn(Float32, nx, ny, n_in, batchsize)
    X1 = randn(Float32, nx, ny, n_in, batchsize)

    # Layer
    H1 = HyperbolicLayer(nx, ny, n_in, batchsize, k, s, p; action=action, α=.2f0, hidden_factor=n_hidden)

    Y0, Y1 = H1.forward(X0, X1)
    X0_, X1_ = H1.inverse(Y0, Y1)

    # Test invertibility
    @test isapprox(norm(X0 - X0_)/norm(X0_), 0f0; atol=1f-6)
    @test isapprox(norm(X1 - X1_)/norm(X1_), 0f0; atol=1f-6)

end

# Loss
function loss(H, X_prev_in, X_curr_in, Y_curr, Y_new)
    X_curr, X_new = H.forward(X_prev_in, X_curr_in)
    f = .5f0*norm(X_new - Y_new)^2 + .5*norm(X_curr - Y_curr)^2
    ΔX_curr = X_curr - Y_curr
    ΔX_new = X_new - Y_new
    ΔX_prev_, ΔX_curr_, X_prev_, X_curr_ = H.backward(ΔX_curr, ΔX_new, X_curr, X_new)

    # Check inverse is correct
    @test isapprox(norm(X_prev_ - X_prev_in)/norm(X_prev_in), 0f0; atol=1f-6)
    @test isapprox(norm(X_curr_ - X_curr_in)/norm(X_curr_in), 0f0; atol=1f-6)

    return f, ΔX_prev_, ΔX_curr_, H.W.grad, H.b.grad
end

function test_grad_X(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    H1 = HyperbolicLayer(nx, ny, n_in, batchsize, k, s, p; action=action, α=.2f0, hidden_factor=n_hidden)
    # Initial guess and residual
    X00 = randn(Float32, nx, ny, n_in, batchsize)
    X01 = randn(Float32, nx, ny, n_in, batchsize)
    dX0 = randn(Float32, nx, ny, n_in, batchsize)
    dX1 = randn(Float32, nx, ny, n_in, batchsize)

    # Observed data
    Y_curr, Y_new = H1.forward(X00, X01)

    # Gradient test for X
    maxiter = 10
    print("\nGradient test hyperbolic layer  input with action=$(action)\n")
    f0, ΔX_prev, ΔX_curr = loss(H1, X00, X01, Y_curr, Y_new)[1:3]
    h = .1f0
    err1 = zeros(Float32, maxiter)
    err2 = zeros(Float32, maxiter)
    for j=1:maxiter
        f = loss(H1, X00 + h*dX0, X01 + h*dX1, Y_curr, Y_new)[1]
        err1[j] = abs(f - f0)
        err2[j] = abs(f - f0 - h*dot(dX0, ΔX_prev) - h*dot(dX1, ΔX_curr))
        print(err1[j], "; ", err2[j], "\n")
        h = h/2f0
    end

    @test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)
end

function test_grad_par(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    # Gradient test for W and b
    H0 = HyperbolicLayer(nx, ny, n_in, batchsize, k, s, p; action="down", α=.2f0, hidden_factor=8)
    Hini = deepcopy(H0)

    dW = randn(Float32, size(H0.W.data))
    db = randn(Float32, size(H0.b.data))
    maxiter = 10

    # Data
    X0 = randn(Float32, nx, ny, n_in, batchsize)
    X1 = randn(Float32, nx, ny, n_in, batchsize)

    # Observed data
    Y_curr, Y_new = H0.forward(X0, X1)

    print("\nGradient test hyperbolic layer weights with action=$(action)\n")
    f0, ΔX_prev, ΔX_curr, ΔW, Δb = loss(H0, X0, X1, Y_curr, Y_new)

    h = .1f0
    err3 = zeros(Float32, maxiter)
    err4 = zeros(Float32, maxiter)

    for j=1:maxiter
        H0.W.data = Hini.W.data + h*dW
        H0.b.data = Hini.b.data + h*db
        f = loss(H0, X0, X1, Y_curr, Y_new)[1]
        err3[j] = abs(f - f0)
        err4[j] = abs(f - f0 - h*dot(dW, ΔW) - h*dot(db, Δb))
        print(err3[j], "; ", err4[j], "\n")
        h = h/2f0
    end
    @test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

end

for action in ["up", "down", "same"]
    test_inv(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    test_grad_X(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    test_grad_par(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
end


text:
# NOTE: None of the models and function in this file are functional. All of these are used for testing purposes and do not represent valid examples either to develop custom
# models. Please refer to the documentation.

struct MockOperationProblem <: PSI.DecisionProblem end
struct MockEmulationProblem <: PSI.EmulationProblem end

function PSI.DecisionModel(
    ::Type{MockOperationProblem},
    ::Type{T},
    sys::PSY.System;
    name=nothing,
    kwargs...,
) where {T <: PM.AbstractPowerModel}
    settings = PSI.Settings(sys; kwargs...)
    return DecisionModel{MockOperationProblem}(
        ProblemTemplate(T),
        sys,
        settings,
        nothing,
        name=name,
    )
end

function make_mock_forecast(horizon, resolution, interval, steps)
    init_time = DateTime("2024-01-01")
    timeseries_data = Dict{Dates.DateTime, Vector{Float64}}()
    for i in 1:steps
        forecast_timestamps = init_time + interval * i
        timeseries_data[forecast_timestamps] = rand(horizon)
    end
    return Deterministic(name="mock_forecast", data=timeseries_data, resolution=resolution)
end

function make_mock_singletimeseries(horizon, resolution)
    init_time = DateTime("2024-01-01")
    tstamps = collect(range(init_time, length=horizon, step=resolution))
    timeseries_data = TimeArray(tstamps, rand(horizon))
    return SingleTimeSeries(name="mock_timeseries", data=timeseries_data)
end

function PSI.DecisionModel(::Type{MockOperationProblem}; name=nothing, kwargs...)
    sys = System(100.0)
    add_component!(sys, Bus(nothing); skip_validation=true)
    l = PowerLoad(nothing)
    gen = ThermalStandard(nothing)
    set_bus!(l, get_component(Bus, sys, "init"))
    set_bus!(gen, get_component(Bus, sys, "init"))
    add_component!(sys, l; skip_validation=true)
    add_component!(sys, gen; skip_validation=true)
    forecast = make_mock_forecast(
        get(kwargs, :horizon, 24),
        get(kwargs, :resolution, Hour(1)),
        get(kwargs, :interval, Hour(1)),
        get(kwargs, :steps, 2),
    )
    add_time_series!(sys, l, forecast)

    settings = PSI.Settings(sys; horizon=get(kwargs, :horizon, 24))
    return DecisionModel{MockOperationProblem}(
        ProblemTemplate(CopperPlatePowerModel),
        sys,
        settings,
        nothing,
        name=name,
    )
end

function PSI.EmulationModel(::Type{MockEmulationProblem}; name=nothing, kwargs...)
    sys = System(100.0)
    add_component!(sys, Bus(nothing); skip_validation=true)
    l = PowerLoad(nothing)
    gen = ThermalStandard(nothing)
    set_bus!(l, get_component(Bus, sys, "init"))
    set_bus!(gen, get_component(Bus, sys, "init"))
    add_component!(sys, l; skip_validation=true)
    add_component!(sys, gen; skip_validation=true)
    single_ts = make_mock_singletimeseries(
        get(kwargs, :horizon, 24),
        get(kwargs, :resolution, Hour(1)),
    )
    add_time_series!(sys, l, single_ts)

    settings = PSI.Settings(sys; horizon=get(kwargs, :horizon, 24))
    return EmulationModel{MockEmulationProblem}(
        ProblemTemplate(CopperPlatePowerModel),
        sys,
        settings,
        nothing,
        name=name,
    )
end

# Only used for testing
function mock_construct_device!(
    problem::PSI.DecisionModel{MockOperationProblem},
    model;
    built_for_recurrent_solves=false,
)
    set_device_model!(problem.template, model)
    template = PSI.get_template(problem)
    PSI.init_optimization_container!(
        PSI.get_optimization_container(problem),
        PSI.get_network_formulation(template),
        PSI.get_system(problem),
    )
    PSI.get_optimization_container(problem).built_for_recurrent_solves =
        built_for_recurrent_solves
    PSI.initialize_system_expressions!(
        PSI.get_optimization_container(problem),
        PSI.get_network_formulation(template),
        PSI.get_system(problem),
    )
    if PSI.validate_available_devices(model, PSI.get_system(problem))
        PSI.construct_device!(
            PSI.get_optimization_container(problem),
            PSI.get_system(problem),
            PSI.ArgumentConstructStage(),
            model,
            PSI.get_network_formulation(template),
        )
        PSI.construct_device!(
            PSI.get_optimization_container(problem),
            PSI.get_system(problem),
            PSI.ModelConstructStage(),
            model,
            PSI.get_network_formulation(template),
        )
    end

    PSI.check_optimization_container(PSI.get_optimization_container(problem))

    JuMP.@objective(
        PSI.get_jump_model(problem),
        MOI.MIN_SENSE,
        PSI.get_objective_fuction(
            PSI.get_optimization_container(problem).objective_function,
        )
    )
    return
end

function mock_construct_network!(problem::PSI.DecisionModel{MockOperationProblem}, model)
    PSI.set_network_model!(problem.template, model)
    PSI.construct_network!(
        PSI.get_optimization_container(problem),
        PSI.get_system(problem),
        model,
        problem.template.branches,
    )
    return
end

function mock_uc_ed_simulation_problems(uc_horizon, ed_horizon)
    return SimulationModels([
        DecisionModel(MockOperationProblem; horizon=uc_horizon, name="UC"),
        DecisionModel(
            MockOperationProblem;
            horizon=ed_horizon,
            resolution=Minute(5),
            name="ED",
        ),
    ])
end

function create_simulation_build_test_problems(
    template_uc=get_template_standard_uc_simulation(),
    template_ed=get_template_nomin_ed_simulation(),
    sys_uc=PSB.build_system(PSITestSystems, "c_sys5_uc"),
    sys_ed=PSB.build_system(PSITestSystems, "c_sys5_ed"),
)
    return SimulationModels(
        decision_models=[
            DecisionModel(template_uc, sys_uc; name="UC", optimizer=GLPK_optimizer),
            DecisionModel(template_ed, sys_ed; name="ED", optimizer=GLPK_optimizer),
        ],
    )
end

struct MockStagesStruct
    stages::Dict{Int, Int}
end

function Base.show(io::IO, struct_stages::MockStagesStruct)
    println(io, "mock problem")
    return
end

function setup_ic_model_container!(model::DecisionModel)
    # This function is only for testing purposes.
    if !PSI.is_empty(model)
        PSI.reset!(model)
    end

    PSI.init_optimization_container!(
        PSI.get_optimization_container(model),
        PSI.get_network_formulation(PSI.get_template(model)),
        PSI.get_system(model),
    )

    PSI.init_model_store_params!(model)

    PSI.populate_aggregated_service_model!(PSI.get_template(model), PSI.get_system(model))
    PSI.populate_contributing_devices!(PSI.get_template(model), PSI.get_system(model))
    PSI.add_services_to_device_model!(PSI.get_template(model))

    @info "Make Initial Conditions Model"
    PSI.set_output_dir!(model, mktempdir(cleanup=true))
    PSI.build_initial_conditions!(model)
    PSI.initialize!(model)
    return
end


text:
# Autogenerated wrapper script for pandoc_jll for x86_64-w64-mingw32
export pandoc

JLLWrappers.@generate_wrapper_header("pandoc")
JLLWrappers.@declare_executable_product(pandoc)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        pandoc,
        "bin\\pandoc.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()