text:
@testset "Unions" begin
    @testset "Low-level" begin
        test(StructDefinition{false}, union_by_name, :VkClearColorValue, :(
            struct _ClearColorValue <: VulkanStruct{false}
                vks::VkClearColorValue
            end
        ))

        test(StructDefinition{false}, union_by_name, :VkClearValue, :(
            struct _ClearValue <: VulkanStruct{false}
                vks::VkClearValue
            end
        ))
    end

    @testset "High-level" begin
        test(StructDefinition{true}, union_by_name, :VkClearColorValue, :(
            struct ClearColorValue <: HighLevelStruct
                vks::VkClearColorValue
            end
        ))

        test(StructDefinition{true}, union_by_name, :VkClearValue, :(
            struct ClearValue <: HighLevelStruct
                vks::VkClearValue
            end
        ))
    end

    @testset "Constructors" begin
        consts = VulkanGen.constructors(StructDefinition{true}(union_by_name(:VkClearValue)))
        expected = [
            :(ClearValue(color::ClearColorValue) = ClearValue(VkClearValue(color.vks))),
            :(ClearValue(depth_stencil::ClearDepthStencilValue) = ClearValue(VkClearValue(_ClearDepthStencilValue(depth_stencil).vks))),
        ]
        foreach(Base.splat(test_ex), zip(consts, expected))
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using Plots, DataFrames, StatsPlots, TimerOutputs

"""
Initialise the basics to avoid massive, incomprehensible globs of code
"""

function fds(func, x, h)
    return (func(x + h) - func(x))/(h)
end

function abs_error(func, x, h, approx_func, accur_value)
    return abs(accur_value - approx_func(func, BigFloat(x), BigFloat(h)))
end

function cds(func, x, h)
    return (func(x + h/2) - func(x - h/2))/(h)
end

function cvm(func, x, h)
    return imag(func(complex(x, h))/h)
end

# *avoid*, not *remove*
sine_set = (x->(sin(x^2)), BigFloat(0.5), big"0.968912421710644784144595449494")
e_set = (x->(ℯ^x), BigFloat(1), ℯ)
atan_set = (x -> (atan(x)/(1 + ℯ^(-x^2))), BigFloat(2), big"0.2746237281548575890153809496775478")

# Test the performance of each approximation for each function

performance_testing_range = map(x->(10.0^(-x)), 3.0:0.01:60)
to = TimerOutput()

optimal_value_df = DataFrame(function_name=String[], approximation_method=String[], optimal_h=BigFloat[], minimal_difference=BigFloat[])

for func_set in [(atan_set, "atan function"), (sine_set, "sine function"), (e_set, "exponent function")]
    for approx_method in [fds, cds, cvm]
        # Simply times how long it takes to collate all absolute errors and find their min
        @timeit to func_set[2]*": "*String(Symbol(approx_method)) begin
            b = map(h->abs_error(func_set[1][1], func_set[1][2], BigFloat(h), approx_method, func_set[1][3]), performance_testing_range)
            c = argmin(b)
        end
        push!(optimal_value_df, (func_set[2], String(Symbol(approx_method)), performance_testing_range[c], b[c]))
    end
end

println(to)
println(optimal_value_df)

# Demonstrate the differences via a chart

charting_range = map(x->(10.0^(-x)), 3.0:0.01:60) # to avoid complete computer death

# Fill the DataFrame with data points for each approximation
df = DataFrame(a=charting_range)
approx_methods_for_df = zip('b':'d', [fds, cds, cvm])
for approx_method ∈ approx_methods_for_df
    df[!, string(approx_method[1])] = map(h->abs_error(atan_set[1], BigFloat(atan_set[2]), BigFloat(h), approx_method[2], atan_set[3]), charting_range)
end

# Plot the DataFrame
inds = (df.c .> 0) .& (df.d .> 0) .& (df.b .> 0)
@df df plot(
    :a[inds], [:b[inds] :c[inds] :d[inds]],
    xaxis=:log10,
    xticks=map(x->(10.0^(-x)), 3.0:10:100.0),
    xflip=true,
    xtickfont = font(5, "Courier"),
    xlabel = "h",
    ylabel = "Absolute Error",
    yscale=:log10,
    label = ["Forward Diff. Scheme" "Central Diff. Scheme" "Complex Var. Method"],
    title = "Demonstration of Approximation Methods"
)

# The differences in time performance for each approximation method were negligible, but the Complex Variable Method appeared to use more significantly more memory. All methods are perfectly capable of approximating a derivative.

# The main difference between these methods is apparent from the graph below. As we can see, the Central Difference Scheme and the Forward Difference Scheme both 'peter out' at sufficiently small values of $h$; after reaching a certain point, the absolute difference for both of these methods begins to rise with a logarithmic trend. This is not the case for the Complex Variable Method - its absolute error continues to plummet logarithmically, then comes to a hard stop within a small range of a minimum. This hard stop remains constant, which allows one to take extremely small values of $n$ without concern for rounding errors.

ext: jl


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ b286fb97-61ed-45d8-9d4c-fea506061e19
using PlutoUI

# ╔═╡ c17b9560-eee8-4216-94fb-fe8db636c03c
md"""

# Introduction to Standards

Standards are at the core of health informatics. However, before diving into the standards of health domain, we are going to start by exploring how standards are an important part of our daily life.

Please  start with [this](https://www.atlasobscura.com/articles/hidden-wonder-international-standards-organization-iso) article, which I hope will give a nice motivation to the idea of standards.

> Standards are like toothbrushes, everyone agrees that they’re a good idea but nobody wants to use anyone else’s. (Attributed to Murtha Baca, Getty Research Institute. Pomerantz, Jeffrey. *Metadata* (MIT Press Essential Knowledge series) (p. 65)).

"""

# ╔═╡ fb66563a-0a13-4dad-8dd8-cc310a6011fa
html"""
Standards Predicted Brexit Decades Ago

"""

# ╔═╡ ba54c574-1507-4528-8542-56f16e8a6bf1
md"""
## What is a standard?

Consulting the Oxford English Dictionary there are, as usual, multiple definitions for _standard_. For our purposes the most relevant definition is

> The authorized exemplar of a unit of measure or weight; e.g. a measuring rod of unit length; a vessel of unit capacity, or a mass of metal of unit weight, preserved in the custody of public officers as a permanent evidence of the legally prescribed magnitude of the unit.

The International Organization for Standards (ISO) provides the following definition of a standard:

> a __document,__ established by consensus and approved by a recognized body, that provides, for common and repeated use, rules, guidelines or characteristics for activities or their results, aimed at the achievement of the optimum degree of order in a given context. ([ISO “What are standards and how do they help?”](https://www.iso.org/sites/ConsumersStandards/1_standards.html#section1_1))

### Key Concept!

The key phrase I want you to keep in mind as we discuss standards is __"aimed at the achievement of the optimum degree of order in a given context."__

As you work through the various sections, think about what the __context__ is and how standards help (or could help) achieve __optimum degree of order__ in that context.

## Standards Everywhere!

My ability to create and your ability to consume these lesson materials are deeply rooted in standards. You do not need to know what kind of computer I'm writing with and I do not need to know what kind of computer you are reading this on, because standards have been created for describing how to save and read characters on a computer. Further you do not need to know how my computer connected to the internet or how the document I created at home was transmitted from my computer to the computer(s) on which Canvas runs.

This ease is not intrinsic to computing but has been created through careful engineering and adoption of standards. Those of you old enough to remember the incompatibility between Mac and PC documents and storage media might have an easier time appreciating this.

English is the [official language of instruction](Your lectures, tutorials, exams, class discussions and other activities will be in English) at the University of Melbourne, but unlike the computing standards described above natural languages like English are more governed by convention. While France has the [Académie Française](https://en.wikipedia.org/wiki/Acad%C3%A9mie_Fran%C3%A7aise) to provide an official, though often still ignored, version of French, there is no equivalent in English.

## Standards are Driven By Technology and Politics

## Standards Examples: Length

Some of the most pervasive standards in our daily life relate to __measurements__.

The English had a unified measuring system data back to the Magna Carta (1215 CE) with constiutant components dating back to Anglo-Saxon and Roman measures. Within [English](https://en.wikipedia.org/wiki/English_units) derived measurements systems, such as [Imperial](https://en.wikipedia.org/wiki/Imperial_units) and [United States Customary Units](https://en.wikipedia.org/wiki/United_States_customary_units), unit bases varied based on application: for example, 12 for length (12 inches in a foot), 16 for volume (16 fluid ounces in a pint), and 16 ounces in a pound. But with lots of other factors joining in: 3 feet in a yard, 14 pounds in a stone,

The French, by contrast, had no common measuring system. After the French Revolution, there was a movement to rationalize French society (as well as to eliminate any connections to the _Ancien Régime_, including measurement, which resulted in the creation of the metric system. In 1795 the meter was [“defined as being one ten millionth of the distance between the North Pole and the Equator through Paris”](https://en.wikipedia.org/wiki/Units_of_measurement_in_France)

As we are probably all aware, the metric system is based on powers of 10. There is nothing special about 10, other than that we have (usually) 10 fingers and (again usually) 10 toes. So there was a _standardization_ of our bases---no more worryiong about whether we need to be thinking in terms of 3, 4, 8, 14, 16, etc. However, from a practical perspective, there are problems with using 10 as a basis---fractions! And absent computers to do teh work for us, humans generally find fractions challenging.

- Few numbers divide evenly into 10 (2 and 5) meaning that you are often having to work with difficult fractions.
	- By contrast, 12 (as in 12 inches in a foot) is evenly divided by 2, 3, 4, and 6.
    - The Babylonians used 60 as their basis. 60 is divided evently by 2, 3, 4, 5, 6, 10, 12, 15, 20, and 30.

Our move from using our fingers for arithmetic to using binary-based computers (computers that count based on powers of 2), introduces new challenges to the metric system since 10 is not a power of 2. So for example,

- One tenth (0.1) cannot be represented exactly in a binary computer, resulting in various computational errors if 0.1 (and similar numbers) form the basis of the computations.

### Moral

Standards are rooted in social systems, solve certain problems (imperfectly), while introducing other problems.

"""

# ╔═╡ f3088148-1cb2-4647-8a9d-71a385aee2da
html"""

You might find it interesting to see who did and who did not sign the original Treaty of the Meter in 1875:

Argentina	 Austria-Hungary	Belgium

Brazil	Denmark	France

Germany 	Italy 	 Peru

Portugal 	 Russia 	 Spain

Sweden and Norway 	 Switzerland 	 Turkey

United States of America 	 Venezuela

"""

# ╔═╡ ec205f4f-c19c-4e99-9f52-f71bc5cbd0df
md"""
## Standards to Standards

As long as relationships between standards are well defined, there is little difference between using one standard or another---the conversions between them is a simple computation.

### Example: US Foot

In the United States, the foot (and hence the inch) has been defined in terms of the meter since 1893. So a foot is __defined as__ 1200/3937 $(1200/3937) meters.

### Example: mmHg

The SI standard for pressure is the Pascal (Pa). However, in medicine blood pressure is reported in millimeters of mercury (mmHG or mm[Hg]), a unit rooted in the historical roots of  the sphygmomanometer.

No worries!

mmHg is __defined as__ exactly 133.322387415 Pa
"""

# ╔═╡ c1a1a8f5-2a75-4eed-b4f4-97bbad2304a0
begin
	solar_year = 365.2422;
	lunar_month = 29.5306;
md"""## Standards Example: Time

Time is a very tricky subject.

- The solar year is approximately $(solar_year) days
- The lunar month is approximately $(lunar_month) days

Divide the solar year by the lunar month and you get
	$(solar_year/lunar_month)

Which is likely where we get our 12 month year, with lots of complicated rules to try to account for the 0.368262073916549 fractional days.

Our 12 month year is likely also the reason we have a 24 hour day, a division created by ancient Babylonian astronomers but not used in common life.

### Technology Determines and Limits Needs for Standards

The standard ancient technology for measuring time was the sundial. The sundial did not lead to our concept of a uniform hour, since the shadow moves slower at mid-day than it does in the evening or the morning. Further, there was probably little social need for small, uniform divisions of time. It was only in the move from an agrarian to an industrial society where shift work needed to be coordinated that a uniform hour started to be defined. The technology of mechanical clocks which started appearing around 1400 CE provided the technological basis for this. Minutes on mechanical clocks came around 200 years later. Seconds arrived on mechanical clocks at nearly the same time.

The classical second was the day dvided by 60x60 (there is that 60 again!). In the 1940s the second was defined as "1/86400" of a mean solar day." ([Wikipedia](https://en.wikipedia.org/wiki/Second#History_of_definition))

## Trains and time tandards

Prior to trains, travel was slow enough that the details of _time-zones_ were not important: it didn't matter whether Philadelphia and New York City, to use a USA example, had the same noon. However, with the speed of trains and the need for passangers to switch trains, coordinating time zones became important for the "optimum degree of order" for travellers.

## Are time standards settled?

Not at all!

In the next section we will learn how seconds are now defined, but since Einstein we have known that the time experienced by two objects (e.g., you and me) is relative to the gravity the objects are experiencing as well as the relative velocity between the objects.

How in the world would this matter? At the very least just think of how dependent we are on GPS satellites, which are experiencing a much weaker gravity than our cars.

### Reflect on
"""
end

# ╔═╡ c08225f3-f69e-4d80-9085-1bc771f341e7
md"""
## [Current Measurement Standards](https://physics.nist.gov/cuu/Units/current.html)

In 2018 our system of measurement standards was completely redesigned based on seven fundamental physical constants.

![NIST SI Constants Cards](https://www.nist.gov/sites/default/files/styles/960_x_960_limit/public/images/2018/10/30/17pml015_si-constants-card.jpg?itok=Ecezgy_j)
NIST SI Constants Cards

- "The meter (m) is defined by taking the fixed numerical value of the speed of light in vacuum c to be 299,792,458 when expressed in the unit m s−1, where the second is defined in terms of ∆νCs."
- "The second, symbol s, is the SI unit of time. It is defined by taking the fixed numerical value of the cesium frequency ΔνCs, the unperturbed ground-state hyperfine transition frequency of the cesium 133 atom, to be 9 192 631 770 when expressed in the unit Hz, which is equal to s-1."

Notice that there is no longer any reference to a day or a year or anything else solar or the equator, Paris or the North Pole.

"""

# ╔═╡ a25033dc-7273-4dd3-a02a-fe3d835523d4
md"""
## The moral from this historical review

> __Standardization becomes critical in a connected world.__ If you never come into contact with your faraway neighbors, it doesn’t especially matter whether you use the same measurements, operate by the same clock, or have the same size windows in your workplace. The Industrial Revolution, however, had the effect of shrinking the world—journeys which would once have been unthinkably long could now take place in the space of weeks, or even days. And for that, there needed to be cooperation along the line.

Similarly, standardization becomes important when the granularity of our interactions with the world increases.

## Standardization in healthcare

Remember the purpose of standards is to achieve "the optimum degree of order in a given context" and that "standardization becomes critical in a connected world."

Our current models of healthcare delivery date from the mid-20th Century. What are some of the ways you think the connectivity and granularity has evolved to make standards more important?

#### Type your answers here
$(@bind medical_connectivity TextField((80,7)))
"""

# ╔═╡ 33f7db73-042b-4b3c-b3f9-e54502a977d5
md"""

## How are standards created?

### There is a standard for that!

Standards are created by __standards developing organizations,__ which are accredited by __national standards organizations__, which join together in __international standards organizations__ to coordinate and published standards of international value.

### SDOs

Standards are created by __Standards Developing Organizations (SDO)__. SDOs can be

- Business groups.
  - Example _USB-DF:_ This SDO consists of businesses that develop USB cables
- Professional organizations.
  - Example _Institute of Electrical and Electroncis Engineers (IEEE):_ Professional organization developing standards for computing, communications, etc.
- Government agencies
  - Example [_Australian Digital Health Agency_](https://www.digitalhealth.gov.au/): Creates standards for Australian-specific needs in the digital health domain.
    - [HL7-FHIR Australia ](https://hl7.com.au/)

### National Standards Organizations (NSO)s

NSOs accredit SDOs within sovereign countries and thus define the rules for standards development. The NSOs are also the representatives for each country in the international standards organizations, such as ISO.

- [Standards Australia](https://www.standards.org.au/)
- [THE AMERICAN NATIONAL STANDARDS INSTITUTE (ANSI)](https://www.ansi.org/)

### International standards organizations

Organizes, publishes.

![standards_bodies.png](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/standards_bodies.png?raw=true)

## Important Medical Standards

- SNOMED-CT

>Today SNOMED International remains a division of CAP and is an American National Standards Institute (ANSI) accredited Standards Developing Organization (SDO), responsible for the continued development and maintenance of the SNOMED Clinical Terms (SNOMED CT) coding standard. ([SNOMED](https://digital.ahrq.gov/key-topics/snomed-ct))

- HL7 (including HL7-FHIR)

>Founded in 1987, Health Level Seven International (HL7) is a not-for-profit, ANSI-accredited standards developing organization dedicated to providing a comprehensive framework and related standards for the exchange, integration, sharing and retrieval of electronic health information that supports clinical practice and the management, delivery and evaluation of health services." ([About HL7 International](https://www.hl7.org/))

- LOINC

>LOINC, which stands for Logical Observation Identifiers Names and Codes, was initiated in 1994 by Clem McDonald, then an investigator at Regenstrief Institute, a non-profit medical research organization associated with Indiana University. Indianapolis-based Regenstrief organized the LOINC Committee to develop a common terminology for laboratory and clinical observations because there was a growing trend to send clinical data electronically from laboratories and other data producers to hospitals, physician's offices, and payers who use the data for clinical care and management purposes.
>
>At the time—and still today—most laboratories and clinical services use HL7 to send their results electronically from their reporting systems to their care systems. However, the tests in these messages are identified by means of their internal, idiosyncratic code values. As a result, receiving care system cannot fully "understand" and properly file the results they receive unless they either adopt the producer's test codes (which is impossible if they receive results from multiple sources), or invest in the work to map each result producer's code system to their internal code system. ([About LOINC](https://loinc.org/about/))

"""

# ╔═╡ 80a39d08-2a1d-428c-bfc4-b505e7050800
md"""
# Health Standards

Digital health initiatives have generally fallen into two categories:

1. How we represent health data/information
   - SNOMED
   - LOINC
1. How we transmit health data/information.
   - HL7

Here we will limit ourselves to data representation, which is the most important concept for our data to knowledge topic.

### Natural language as representation

Much of human language addresses representing data in the world as well as the abstract  concepts we create to provide meaning about the data. Healthcare makes extensive use of natural language to represent and communicate data. However, natural language is a far from perfect medium for data representation. Here are some problems with natural langague:

-  Polysemy (homograph): A single term refers to multiple concepts. The word _set_ has [430 senses listed in the Oxford English Dictionary](https://www.guinnessworldrecords.com/world-records/english-word-with-the-most-meanings)
   - Mercury could refer to a mythical being, a planet, a metal, an automobile
-  Synonymy: Multiple terms refer to the same concept
   – AMI=Acute Myocardial Infarctus=Infarctus of the myocardium=Heart attack
   - salinity=saltiness
-  Imprecision: Lack of specificities. Natural speech often is elliptical relying on the knowledge of the listener to fill in the missing details. This can be problematic.
   – e.g., infarctus?

These all lead to __ambiguity__, which is not conducive to a "optimum degree of order."

## Standards Solution

Most of the standards that we will examine here are [ontologies](https://en.wikipedia.org/wiki/Ontology_(information_science)). Ontologies model the concepts of a domain using a heirarchy of concepts. That is, concepts have "is a" and "has a" relationships. As an illustration of this, the figure below shows a standard (Uber ontology) representation of the right kidney.

![dimmed kidney](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/uber_kidney_dimmed.png?raw=true)

We can see that the _right kidney_ is a kind of _kidney_ which is a kind of _cavitated compound organ_ all the way up to an _entity_. In this figure I've dimmed some of the nodes in the representation so that only a single pathway is shown. The next figure shows the complete representation.

![uber full kidney](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/uber_kidney.png?raw=true)

Here you can see that a kidney can be viewed in a variety of ways:

1. It is a an organ
1. It is an element of the abdomen
1. It is a lateral anatomical structure.

This is an example of a poly-heirarchy, which we will touch on again below.

### Concept-Base Representation

The foundational idea for standards-based solutions is referred to as "Concept Based Representation" (CBR).

![Concept based representation](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/concept_based_representation.png?raw=true)

CBR has three parts

1. A unit of thought (concept)
1. A unique identifier for that unit of thought
1. A lexical expression (term) for that unit of thought

### A Few Principles of Successful Standards

1. Formal definition
   - The standard provides a formal definition of the concept (unit of thought) being represented.
1. Terminology and synonyms
   - The standard provides the preferred term for the concept as well as any "allowed" synonyms
1. Concept orientation
   - The standard ensures that each concept has one unambiguous, non-redundant meaning
1. Concept permanence
   - Old concepts cannot be deleted from the standard.
1. Non-semantic identifiers
   - The standard provides an ID for the concept that is not a name, that is not a code that will run out of room, and that is not  simplya place in a hierarchy.
1. Poly-hierarchy
   - The standard allows for a concept to have multiple parent concepts.

### Example: Lisinopril

We will illustrate some of these concepts using the drug _lisinopril_. We will look at lisinopril in the following standards

- [SNOMED-CT](http://purl.bioontology.org/ontology/SNOMEDCT/386873009)
- [RxNORM](http://purl.bioontology.org/ontology/RXNORM/29046)
- [Chemical Entities of Biological Interest Ontology](http://purl.obolibrary.org/obo/CHEBI_43755)
- [MeSH](http://purl.bioontology.org/ontology/MESH/D017706)

- Australian Medicines Terminology
- Australian PBS

#### 1. Formal Definition of `lisinopril`

Only MeSH and Chemical Entities provide a definition of lisinopril.

- Mesh

>One of the ANGIOTENSIN-CONVERTING ENZYME INHIBITORS (ACE inhibitors), orally active, that has been used in the treatment of hypertension and congestive heart failure.

- Chemical Entities

>C21H31N3O5

#### 2. Terminology and synonyms

All the standards provide a preferred name and synonyms. In the case of RxNORM synonyms are "trade names".

#### 5. Non-semantic identifiers

All the standards provide unique, non-semantic identifiers

#### 6. Poly-hierarchy

Only SNOMED-CT has a poly-hierarchy representation of lisinopril. Viewing the graph below shows that SNOMED represents lisinopril both in terms of its chemical composition (dipeptide) as well as its clinical purpose ("hypotensive agent").

![lisinopril SNOMED Visualization](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/snomed_lisinopril_viz.png?raw=true)

The figure is an illustration of how ontologies can help with data analysis. We can see that lisinopril is one of many ACE inhibitors which is a kind of "Renin-angiotensin system inhibitor" which is a kind of "Hypotensive agent." Thus data could be aggregated at various levels of granualrity.

![ACE inhibitors SNOMED](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/snomed_ace_inhibitors.png?raw=true)

### Australian Resources

Because standards are intended for a specific context, Australia has its own set of health standards.

#### SNOMED-CT-AU

Some of these standards are expansion of existing international standards. Australia has its own modifications of SNOMED-CT which can be explored using CSIRO's [Shrimp Browser](https://ontoserver.csiro.au/shrimp).

#### Australian Medicine Terminology

There is also the Australian Medicine Terminology (AMT), something akin to RxNORM that we explored above using Bioportal. The Australian Medicine Terminology is also accessible via Shrimp, however, it is a little hard to find. Below shows the AMT heirarchy of lisinopril 10 mg.

![Shrimp lisinopril 10 mg](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/amt_lisinopril_heirarchy.png?raw=true).

The next figure shows the details of how a particular lisinopril 10 mg concept is represented in ATM

![lisinopril AU](https://github.com/chapmanbe/raw-stuff/blob/main/Standards/lisinopril_aumed1.png?raw=true)

#### Prescription Benefit Scheme

## Why are there so many standards?

Recall "Standards are like toothbrushes, everyone agrees that they’re a good idea but nobody wants to use anyone else’s."

But also standards are developed for a particular context and are abstractions for that context. Recalling Picasso's bulls, the _essence_ of lisinopril differs, depending on the context.

"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─c17b9560-eee8-4216-94fb-fe8db636c03c
# ╟─fb66563a-0a13-4dad-8dd8-cc310a6011fa
# ╠═ba54c574-1507-4528-8542-56f16e8a6bf1
# ╟─f3088148-1cb2-4647-8a9d-71a385aee2da
# ╠═ec205f4f-c19c-4e99-9f52-f71bc5cbd0df
# ╠═c1a1a8f5-2a75-4eed-b4f4-97bbad2304a0
# ╟─c08225f3-f69e-4d80-9085-1bc771f341e7
# ╟─a25033dc-7273-4dd3-a02a-fe3d835523d4
# ╠═33f7db73-042b-4b3c-b3f9-e54502a977d5
# ╠═80a39d08-2a1d-428c-bfc4-b505e7050800
# ╟─b286fb97-61ed-45d8-9d4c-fea506061e19
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
Formulation described in:

    Pan, K., & Guan, Y. (2016). Strong formulations for multistage stochastic
    self-scheduling unit commitment. Operations Research, 64(6), 1482-1498.
    DOI: https://doi.org/10.1287/opre.2016.1520
"""
module PanGua2016

import ..RampingFormulation

struct Ramping <: RampingFormulation end

end


text:
@testset "Limited" begin
    @testset "Limited Array" begin
        L1 = Limited([1, 2, 3]; limit = 5)

        @test n_dropped(L1) == 0

        for i in 4:8
            push!(L1, i)
        end

        @test n_dropped(L1) == 3

        L1[1] = 0

        @test L1[1] == 0
    end

    @testset "Limited Dict" begin
        L1 = Limited(Dict("a" => 1, "b" => 2); limit = 1)

        @test n_dropped(L1) == 1

        L1["c"] = 3

        @test_throws KeyError L1["c"]
        @test n_dropped(L1) == 2
    end
end

@testset "StaticAttrs" begin
    A = StaticAttrs()
    @test length(A) == 0

    A1 = StaticAttrs((x = 1, y = -2.0, z = "abcdefg"); value_length_limit = 3)
    @test A1["x"] == 1
    @test A1["z"] == "abc"

    A2 = StaticAttrs((; :x => 1, :y => -2.0, :z => "0"))
    @test A2[:z] == "0"

    @test_throws MethodError StaticAttrs("abc" => nothing)
end

@testset "DynamicAttrs" begin
    d = DynamicAttrs("a" => 1; count_limit = 4, value_length_limit = 5)
    d["b"] = "bbbbb"
    d["c"] = "ccccc" * "extra"
    d["d"] = ["ddddd" * "extra", "dd"]
    d["e"] = "???"

    @test length(d) == 4
    @test n_dropped(d) == 1
    @test d["b"] == "bbbbb"
    @test d["c"] == "ccccc"
    @test d["d"] == ["ddddd", "dd"]

    d["b"] = "bbbbbbbbbb"
    @test d["b"] == "bbbbb"

    # count_limit exceeded
    d["x"] = "abcdefg"
    @test_throws KeyError d["x"]
end

@testset "context" begin
    @test isnothing(get(current_context(), :say, nothing))

    with_context(; say = "foo") do
        @test current_context()[:say] == "foo"
        with_context(; say = "bar") do
            @test current_context()[:say] == "bar"
        end
        @test current_context()[:say] == "foo"
    end

    @test isnothing(get(current_context(), :say, nothing))

    with_context(; a = "foo") do
        Threads.@spawn begin
            with_context(; b = "bar") do
                Threads.@spawn begin
                    @test current_context()[:a] == "foo"
                    @test current_context()[:b] == "bar"
                end
            end
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, BinaryBuilderBase, Pkg

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.

name = "GAP"
version = v"400.1100.1"
upstream_version = v"4.11.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/gap-system/gap.git", "42b0b58709ccb7272de7ad467ecd6f05f570c817"),
#    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(upstream_version)/gap-$(upstream_version)-core.tar.bz2",
#                  "6637f66409bc91af21eaa38368153270b71b13b55b75cc1550ed867c629901d1"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gap*

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/configure.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/Makefile.patch

# run autogen.sh if compiling from it source and/or if configure was patched
./autogen.sh

# compile GAP
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-gmp=${prefix} \
    --with-readline=${prefix} \
    --with-zlib=${prefix} \
    --with-gc=julia \
    --with-julia
make -j${nproc}

# install GAP binaries
make install-bin install-headers install-libgap

# also install config.h
cp gen/config.h ${prefix}/include/gap

# get rid of *.la files, they just cause trouble
rm ${prefix}/lib/*.la

# get rid of the wrapper shell script, which is useless for us
mv ${prefix}/bin/gap.real ${prefix}/bin/gap

# install gac and sysinfo.gap
mkdir -p ${prefix}/share/gap/
cp gac sysinfo.gap ${prefix}/share/gap/

# We deliberately do NOT install the GAP library, documentation, etc. because
# they are identical across all platforms; instead, we use another platform
# independent artifact to ship them to the user.
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# we only care about 64bit builds
filter!(p -> nbits(p) == 64, platforms)

# Windows is not supported
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("gap", :gap),
    LibraryProduct("libgap", :libgap),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("Readline_jll"),
    Dependency("Zlib_jll"),

    # GAP tries hard to produce a binary that works in all Julia versions,
    # regardless of which version of Julia it was compiled again; so the
    # version restriction below could be dropped or changed if necessary
    BuildDependency(PackageSpec(name="libjulia_jll", version=v"1.5.3")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"7",
               init_block = """

    sym = dlsym(libgap_handle, :GAP_InitJuliaMemoryInterface)
    ccall(sym, Nothing, (Any, Ptr{Nothing}), @__MODULE__, C_NULL)
""")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# adapted from http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw4b.py

import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    # version 1
    @time GU.read("ospy/data4/aster.img") do ds
        count = 0
        total = 0
        data = GU.fetch(ds, 1)
        for (cols,rows) in GU.WindowIterator(GU.fetchband(ds, 1))
            window = data[cols, rows]
            count = count + sum(window .> 0)
            total = total + sum(window)
        end
        println("Ignoring 0:  $(total / count)")
        println("Including 0: $(total / (GU.height(ds) * GU.width(ds)))")
    end
    # Ignoring 0: 76.33891347095299
    # Including 0: 47.55674749653172
    #   0.367471 seconds (173.86 k allocations: 91.463 MB, 3.39% gc time)

    # version 2
    @time GU.read("ospy/data4/aster.img") do ds
        band = GU.fetchband(ds, 1)
        count = 0
        total = 0
        buffer = Array(GU.getdatatype(band), GU.getblocksize(band)...)
        for (cols,rows) in GU.WindowIterator(band)
            GU.rasterio!(band, buffer, rows, cols)
            data = buffer[1:length(cols),1:length(rows)]
            count += sum(data .> 0)
            total += sum(data)
        end
        println("Ignoring 0:  $(total / count)")
        println("Including 0: $(total / (GU.height(ds) * GU.width(ds)))")
    end
    # Ignoring 0:  76.33891347095299
    # Including 0: 47.55674749653172
    #   0.337447 seconds (208.89 k allocations: 65.353 MB, 7.78% gc time)

    # version 3
    @time GU.read("ospy/data4/aster.img") do ds
        count = 0
        total = 0
        # BufferIterator uses a single buffer, so this loop cannot be parallelized
        for data in GU.BufferIterator(GU.fetchband(ds, 1))
            count += sum(data .> 0)
            total += sum(data)
        end
        println("Ignoring 0:  $(total / count)")
        println("Including 0: $(total / (GU.height(ds) * GU.width(ds)))")
    end
    # Ignoring 0:  76.33891347095299
    # Including 0: 47.55674749653172
    #   0.312078 seconds (205.92 k allocations: 65.093 MB, 8.70% gc time)

end


text:
## Welfare analysis. INCOMPLETE, MISSING THE 3D GRAPHICS [afaik, the best way to make it in Julia is through Generic Mapping Tools (using GMT.jl)

q=(v0_transition = rt[2].v,
q.v0_transition_t3 = rt[3].v,
q.v0_ss = rt[1].v,
q.v0_ss_final = rt[s.N].v,

#Individual-specific gains
ind_welfare = (q.v0_transition./q.v0_ss).^(1/(1-m.γ)) .- 1,
ind_welfare_t3 = (q.v0_transition_t3./q.v0_ss).^(1/(1-m.γ)) .- 1,
measure_0 = sim_0.measure,
measure_end = sim_end.measure,
ind_welfare_avg = sum(q.measure_0.*q.ind_welfare),

#Gains without transition
ind_welfare_notrans = (q.v0_ss_final./q.v0_ss).^(1/(1-m.γ)) .- 1,
ind_welfare_avg_notrans = sum(q.measure_0.*q.ind_welfare_notrans),

#Average gains
agg_welfare_avg = (sum(q.measure_0.*q.v0_transition)/sum(q.measure_0.*q.v0_ss)).^(1/(1-m.γ)) .- 1,

#Gains without transition
agg_welfare_avg_notrans = (sum(q.measure_end.*q.v0_ss_final)/sum(q.measure_0.*q.v0_ss)).^(1/(1-m.γ)) .- 1,
agg_welfare_avg_notrans_fixedmeasure = (sum(q.measure_0.*q.v0_ss_final)/sum(q.measure_0.*q.v0_ss)).^(1/(1-m.γ)) .- 1,
agg_welfare_avg_notrans_fixedvalue = (sum(q.measure_end.*q.v0_ss)/sum(q.measure_0.*q.v0_ss)).^(1/(1-m.γ)) .- 1)

### STATISTICS

# All
q=merge(q,(share_winners = sum(sim_0.measure[q.ind_welfare.>0]),
share_losers = sum(sim_0.measure[q.ind_welfare.<0]),
welfare_avg_winners = sum(sim_0.measure[q.ind_welfare.>=0].*q.ind_welfare[q.ind_welfare.>=0])/ sum(sim_0.measure[q.ind_welfare.>=0]),
welfare_avg_losers = sum(sim_0.measure[q.ind_welfare.<0].*q.ind_welfare[q.ind_welfare.<0])/ sum(sim_0.measure[q.ind_welfare.<0])))

# Exporters vs Non-exporters
q=merge(q,(welfare_avg_exp = sum(sim_0.measure[rt[1].e.==1].*q.ind_welfare[rt[1].e.==1])/sum(sim_0.measure[rt[1].e.==1]),
welfare_avg_nonexp = sum(sim_0.measure[rt[1].e.==0].*q.ind_welfare[rt[1].e.==0])/sum(sim_0.measure[rt[1].e.==0]),

share_losers_exp = sum(sim_0.measure[q.ind_welfare.<0 .& rt[1].e.==1])/ sum(sim_0.measure[rt[1].e.==1]),
share_losers_nonexp = sum(sim_0.measure[q.ind_welfare.<0 .& rt[1].e.==0])/  sum(sim_0.measure[rt[1].e.==0]),
share_winners_exp = sum(sim_0.measure[q.ind_welfare.>=0 .& rt[1].e.==1]) / sum(sim_0.measure[rt[1].e.==1]),
share_winners_nonexp = sum(sim_0.measure[q.ind_welfare.>=0 .& rt[1].e.==0]) / sum(sim_0.measure[rt[1].e.==0]),

welfare_avg_losers_exp=  sum(sim_0.measure[rt[1].e.==1].*q.ind_welfare[rt[1].e.==1].*(q.ind_welfare[rt[1].e.==1].<0))/sum(sim_0.measure[rt[1].e.==1].*(q.ind_welfare[rt[1].e.==1].<0)),
welfare_avg_winners_exp=  sum(sim_0.measure[rt[1].e.==1].*q.ind_welfare[rt[1].e.==1].*(q.ind_welfare[rt[1].e.==1].>=0))/sum(sim_0.measure[rt[1].e.==1].*(q.ind_welfare[rt[1].e.==1].>=0)),
welfare_avg_losers_nonexp=  sum(sim_0.measure[rt[1].e.==0].*q.ind_welfare[rt[1].e.==0].*(q.ind_welfare[rt[1].e.==0].<0))/sum(sim_0.measure[rt[1].e.==0].*(q.ind_welfare[rt[1].e.==0].<0)),
welfare_avg_winners_nonexp=  sum(sim_0.measure[rt[1].e.==0].*q.ind_welfare[rt[1].e.==0].*(q.ind_welfare[rt[1].e.==0].>=0))/sum(sim_0.measure[rt[1].e.==0].*(q.ind_welfare[rt[1].e.==0].>=0))))

# Rich vs Poor
q=merge(q,(measure_a_cum = cumsum(sum(sim_0.measure,dims=2),dims=1),
aP = sum(q.measure_a_cum.<0.99,dims=1) .+1,
# aP = s.a_grid_size./2 .+1; # threshold for asset

welfare_avg_poor = sum(sim_0.measure[1:q.aP-1,:].*q.ind_welfare[1:q.aP-1,:])/sum(sim_0.measure[1:q.aP-1,:]),
welfare_avg_rich = sum(sim_0.measure[q.aP:end,:].*q.ind_welfare[q.aP:end,:])/sum(sim_0.measure[q.aP:end,:]),

share_losers_poor = sum(sim_0.measure[1:q.aP-1,:].*(q.ind_welfare[1:q.aP-1,:].<0))/sum(sim_0.measure[1:q.aP-1,:]),
share_winners_poor = sum(sim_0.measure[1:q.aP-1,:].*(q.ind_welfare[1:q.aP-1,:].>=0))/sum(sim_0.measure[1:q.aP-1,:]),
welfare_avg_losers_poor =  sum(sim_0.measure[1:q.aP-1,:].*q.ind_welfare[1:q.aP-1,:].*(q.ind_welfare[1:q.aP-1,:].<0))/sum(sim_0.measure[1:q.aP-1,:].*(q.ind_welfare[1:q.aP-1,:].<0)),
welfare_avg_winners_poor =  sum(sim_0.measure[1:q.aP-1,:].*q.ind_welfare[1:q.aP-1,:].*(q.ind_welfare[1:q.aP-1,:].>=0))/sum(sim_0.measure[1:q.aP-1,:].*(q.ind_welfare[1:q.aP-1,:].>=0)),

share_losers_rich = sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].<0))/sum(sim_0.measure[q.aP:end,:]),
share_winners_rich = sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].>=0))/sum(sim_0.measure[q.aP:end,:]),
welfare_avg_losers_rich = sum(sim_0.measure[q.aP:end,:].*q.ind_welfare[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].<0))/sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].<0)),
welfare_avg_winners_rich = sum(sim_0.measure[q.aP:end,:].*q.ind_welfare[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].>=0))/sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].>=0))))

# Productive vs Non-productive
q=merge(q,(measure_z_cum = cumsum(sum(sim_0.measure,dims=1),dims=2),
zP = sum(q.measure_z_cum.<0.5,dims=1).+1,
# zP = s.z_grid_size./2 .+1; # threshold for productive

welfare_avg_NP = sum(sim_0.measure[:,1:q.zP-1].*q.ind_welfare[:,1:q.zP-1])/sum(sim_0.measure[:,1:q.zP-1]),
welfare_avg_P = sum(sim_0.measure[:,q.zP:end].*q.ind_welfare[:,q.zP:end])/sum(sim_0.measure[:,q.zP:end]),
share_losers_NP =  sum(sim_0.measure[:,1:q.zP-1].*(q.ind_welfare[:,1:q.zP-1].<0))/sum(sim_0.measure[:,1:q.zP-1]),
share_winners_NP =  sum(sim_0.measure[:,1:q.zP-1].*(q.ind_welfare[:,1:q.zP-1].>=0))/sum(sim_0.measure[:,1:q.zP-1]),
welfare_avg_winners_NP  = sum(sim_0.measure[:,1:q.zP-1].*q.ind_welfare[:,1:q.zP-1].*(q.ind_welfare[:,1:q.zP-1].>=0))/sum(sim_0.measure[:,1:q.zP-1].*(q.ind_welfare[:,1:q.zP-1].>=0)),
welfare_avg_losers_NP  = sum(sim_0.measure[:,1:q.zP-1].*q.ind_welfare[:,1:q.zP-1].*(q.ind_welfare[:,1:q.zP-1].<0))/sum(sim_0.measure[:,1:q.zP-1].*(q.ind_welfare[:,1:q.zP-1].<0)),
share_winners_P=  sum(sim_0.measure[:,q.zP:end].*(q.ind_welfare[:,q.zP:end].>=0))/sum(sim_0.measure[:,q.zP:end]),
share_losers_P =  sum(sim_0.measure[:,q.zP:end].*(q.ind_welfare[:,q.zP:end].<0))/sum(sim_0.measure[:,q.zP:end]),
welfare_avg_losers_P = sum(sim_0.measure[:,q.zP:end].*q.ind_welfare[:,q.zP:end].*(q.ind_welfare[:,q.zP:end].<0))/sum(sim_0.measure[:,q.zP:end].*(q.ind_welfare[:,q.zP:end].<0)),
welfare_avg_winners_P = sum(sim_0.measure[:,q.zP:end].*q.ind_welfare[:,q.zP:end].*(q.ind_welfare[:,q.zP:end].>=0))/sum(sim_0.measure[:,q.zP:end].*(q.ind_welfare[:,q.zP:end].>=0))))

###  Exporter vs Nonexporters and Rich vs Poor
q=merge(q,(measure_a_cum = cumsum(sum(sim_0.measure,dims=2),dims=1),
aP = sum(q.measure_a_cum .<0.99,dims=1) .+1,
# aP = s.a_grid_size/2+1, # threshold for asset
aVP = 1,

welfare_avg_expR =  sum((rt[1].e[q.aP:end,:].==1).*sim_0.measure[q.aP:end,:].*q.ind_welfare[q.aP:end,:])/sum((rt[1].e[q.aP:end,:].==1).*sim_0.measure[q.aP:end,:]),
welfare_avg_nonexpR =  sum((rt[1].e[q.aP:end,:].==0).*sim_0.measure[q.aP:end,:].*q.ind_welfare[q.aP:end,:])/sum((rt[1].e[q.aP:end,:].==0).*sim_0.measure[q.aP:end,:]),
welfare_avg_expP =  sum((rt[1].e[q.aVP+1:q.aP-1,:].==1).*sim_0.measure[q.aVP+1:q.aP-1,:].*q.ind_welfare[q.aVP+1:q.aP-1,:])/sum((rt[1].e[q.aVP+1:q.aP-1,:].==1).*sim_0.measure[q.aVP+1:q.aP-1,:]),
welfare_avg_nonexpP =  sum((rt[1].e[q.aVP+1:q.aP-1,:].==0).*sim_0.measure[q.aVP+1:q.aP-1,:].*q.ind_welfare[q.aVP+1:q.aP-1,:])/sum((rt[1].e[q.aVP+1:q.aP-1,:].==0).*sim_0.measure[q.aVP+1:q.aP-1,:]),
welfare_avg_expVP = sum((rt[1].e[1:q.aVP ,:].==1).*sim_0.measure[1:q.aVP ,:].*q.ind_welfare[1:q.aVP ,:])/sum((rt[1].e[1:q.aVP ,:].==1).*sim_0.measure[1:q.aVP ,:]),
welfare_avg_nonexpVP = sum((rt[1].e[1:q.aVP ,:].==0).*sim_0.measure[1:q.aVP ,:].*q.ind_welfare[1:q.aVP ,:])/sum((rt[1].e[1:q.aVP ,:].==0).*sim_0.measure[1:q.aVP ,:]),

share_expR = sum((rt[1].e[q.aP:end,:].==1).*sim_0.measure[q.aP:end,:]),
share_expP = sum((rt[1].e[q.aVP+1:q.aP-1,:].==1).*sim_0.measure[q.aVP+1:q.aP-1,:]),
share_nonexpR = sum((rt[1].e[q.aP:end,:].==0).*sim_0.measure[q.aP:end,:]),
share_nonexpP = sum((rt[1].e[q.aVP+1:q.aP-1,:].==0).*sim_0.measure[q.aVP+1:q.aP-1,:]),
share_expVP = sum((rt[1].e[1:q.aVP,:].==1).*sim_0.measure[1:q.aVP,:]),
share_nonexpVP = sum((rt[1].e[1:q.aVP,:].==0).*sim_0.measure[1:q.aVP,:]),

share_losers_expR = sum((q.ind_welfare[q.aP:end,:].<0) .*  (rt[1].e[q.aP:end,:].==1).*sim_0.measure[q.aP:end,:]),
share_winners_expR = sum((q.ind_welfare[q.aP:end,:].>=0) .*  (rt[1].e[q.aP:end,:].==1).*sim_0.measure[q.aP:end,:]),
share_losers_nonexpR = sum( (q.ind_welfare[q.aP:end,:].<0) .*  (rt[1].e[q.aP:end,:].==0).*sim_0.measure[q.aP:end,:]),
share_winners_nonexpR = sum( (q.ind_welfare[q.aP:end,:].>=0) .*  (rt[1].e[q.aP:end,:].==0).*sim_0.measure[q.aP:end,:]),

share_losers_expP = sum( (q.ind_welfare[q.aVP+1:q.aP-1,:].<0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==1).*sim_0.measure[q.aVP+1:q.aP-1,:]),
share_winners_expP = sum( (q.ind_welfare[q.aVP+1:q.aP-1,:].>=0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==1).*sim_0.measure[q.aVP+1:q.aP-1,:]),
share_losers_nonexpP = sum( (q.ind_welfare[q.aVP+1:q.aP-1,:].<0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==0).*sim_0.measure[q.aVP+1:q.aP-1,:]),
share_winners_nonexpP = sum((q.ind_welfare[q.aVP+1:q.aP-1,:].>=0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==0).*sim_0.measure[q.aVP+1:q.aP-1,:]),

share_losers_expVP = sum((q.ind_welfare[1:q.aVP,:].<0) .*  (rt[1].e[1:q.aVP,:].==1).*sim_0.measure[1:q.aVP,:]),
share_winners_expVP = sum((q.ind_welfare[1:q.aVP,:].>=0) .*  (rt[1].e[1:q.aVP,:].==1).*sim_0.measure[1:q.aVP,:]),
share_losers_nonexpVP = sum((q.ind_welfare[1:q.aVP,:].<0) .*  (rt[1].e[1:q.aVP,:].==0).*sim_0.measure[1:q.aVP,:]),
share_winners_nonexpVP = sum((q.ind_welfare[1:q.aVP,:].>=0) .*  (rt[1].e[1:q.aVP,:].==0).*sim_0.measure[1:q.aVP,:]),

welfare_avg_losers_expR = sum(sim_0.measure[q.aP:end,:] .* q.ind_welfare[q.aP:end,:] .* (q.ind_welfare[q.aP:end,:].<0) .*(rt[1].e[q.aP:end,:].==1))sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].<0) .*  (rt[1].e[q.aP:end,:].==1)),
welfare_avg_winners_expR = sum(sim_0.measure[q.aP:end,:].*q.ind_welfare[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].>=0) .* (rt[1].e[q.aP:end,:].==1))/sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].>=0) .* (rt[1].e[q.aP:end,:].==1)),
welfare_avg_losers_nonexpR =sum(sim_0.measure[q.aP:end,:] .* q.ind_welfare[q.aP:end,:] .* (q.ind_welfare[q.aP:end,:].<0) .*  (rt[1].e[q.aP:end,:].==0))/sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].<0) .*  (rt[1].e[q.aP:end,:].==0)),
welfare_avg_winners_nonexpR = sum(sim_0.measure[q.aP:end,:].*q.ind_welfare[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].>=0) .* (rt[1].e[q.aP:end,:].==0))/sum(sim_0.measure[q.aP:end,:].*(q.ind_welfare[q.aP:end,:].>=0).*(rt[1].e[q.aP:end,:].==0)),

welfare_avg_losers_expP = sum(sim_0.measure[q.aVP+1:q.aP-1,:] .*q.ind_welfare[q.aVP+1:q.aP-1,:].* (q.ind_welfare[q.aVP+1:q.aP-1,:].<0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==1))/sum(sim_0.measure[q.aVP+1:q.aP-1,:].*(q.ind_welfare[q.aVP+1:q.aP-1,:].<0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==1) ),
welfare_avg_winners_expP = sum(sim_0.measure[q.aVP+1:q.aP-1,:].*q.ind_welfare[q.aVP+1:q.aP-1,:].*(q.ind_welfare[q.aVP+1:q.aP-1,:].>=0) .* (rt[1].e[q.aVP+1:q.aP-1,:].==1))/sum(sim_0.measure[q.aVP+1:q.aP-1,:].*(q.ind_welfare[q.aVP+1:q.aP-1,:].>=0) .* (rt[1].e[q.aVP+1:q.aP-1,:].==1) ),
welfare_avg_losers_nonexpP = sum(sim_0.measure[q.aVP+1:q.aP-1,:] .* q.ind_welfare[q.aVP+1:q.aP-1,:] .* (q.ind_welfare[q.aVP+1:q.aP-1,:].<0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==0))/sum(sim_0.measure[q.aVP+1:q.aP-1,:].*(q.ind_welfare[q.aVP+1:q.aP-1,:].<0) .*  (rt[1].e[q.aVP+1:q.aP-1,:].==0) ),
welfare_avg_winners_nonexpP = sum(sim_0.measure[q.aVP+1:q.aP-1,:].*q.ind_welfare[q.aVP+1:q.aP-1,:].*(q.ind_welfare[q.aVP+1:q.aP-1,:].>=0) .* (rt[1].e[q.aVP+1:q.aP-1,:].==0))/sum(sim_0.measure[q.aVP+1:q.aP-1,:].*(q.ind_welfare[q.aVP+1:q.aP-1,:].>=0) .* (rt[1].e[q.aVP+1:q.aP-1,:].==0) ),

welfare_avg_losers_expVP = sum(sim_0.measure[1:q.aVP,:] .* q.ind_welfare[1:q.aVP,:] .* (q.ind_welfare[1:q.aVP,:].<0) .*  (rt[1].e[1:q.aVP,:].==1))/sum(sim_0.measure[1:q.aVP,:].*(q.ind_welfare[1:q.aVP,:].<0) .*  (rt[1].e[1:q.aVP,:].==1)),
q.welfare_avg_winners_expVP = sum(sim_0.measure[1:q.aVP,:].*q.ind_welfare[1:q.aVP,:].*(q.ind_welfare[1:q.aVP,:].>=0) .* (rt[1].e[1:q.aVP,:].==1))/sum(sim_0.measure[1:q.aVP,:].*(q.ind_welfare[1:q.aVP,:].>=0) .* (rt[1].e[1:q.aVP,:].==1)),
q.welfare_avg_losers_nonexpVP =sum(sim_0.measure[1:q.aVP,:] .* q.ind_welfare[1:q.aVP,:] .* (q.ind_welfare[1:q.aVP,:].<0) .*  (rt[1].e[1:q.aVP,:].==0))/sum(sim_0.measure[1:q.aVP,:].*(q.ind_welfare[1:q.aVP,:].<0) .*  (rt[1].e[1:q.aVP,:].==0)),
q.welfare_avg_winners_nonexpVP = sum(sim_0.measure[1:q.aVP,:].*q.ind_welfare[1:q.aVP,:].*(q.ind_welfare[1:q.aVP,:].>=0) .* (rt[1].e[1:q.aVP,:].==0))/sum(sim_0.measure[1:q.aVP,:].*(q.ind_welfare[1:q.aVP,:].>=0) .* (rt[1].e[1:q.aVP,:].==0))))

###  Productive vs Non-productive and Rich vs Poor
q=merge(q,(zP_low = sum(q.measure_z_cum .<0.33,dims=1) .+1,
zP_med = sum(q.measure_z_cum.<0.66,dims=1) .+1,

# shares
share_lowR = sum(sim_0.measure[q.aP:end,1:q.zP_low.-1]),
share_lowP = sum(sim_0.measure[q.aVP.+1:q.aP-1,1:q.zP_low.-1]),
share_lowVP = sum(sim_0.measure[1:q.aVP,1:q.zP_low.-1]),

share_medR = sum(sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1]),
share_medP = sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1]),
share_medVP = sum(sim_0.measure[1:q.aVP,q.zP_low:q.zP_med]),

share_highR = sum(sim_0.measure[q.aP:end,q.zP_med:end]),
share_highP = sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end]),
share_highVP = sum(sim_0.measure[1:q.aVP,q.zP_med:end])))

# welfare
q=merge(q,(welfare_avg_lowR =  sum(sim_0.measure[q.aP:end,1:q.zP_low.-1].*q.ind_welfare[q.aP:end,1:q.zP_low.-1])/sum(sim_0.measure[q.aP:end,1:q.zP_low.-1]),
welfare_avg_lowP =  sum(sim_0.measure[q.aVP.+1:q.aP.-1,1:q.zP_low.-1].*q.ind_welfare[q.aVP.+1:q.aP.-1,1:q.zP_low.-1])/sum(sim_0.measure[q.aVP.+1:q.aP.-1,1:q.zP_low.-1]),
welfare_avg_lowVP =  sum(sim_0.measure[1:q.aVP ,1:q.zP_low.-1].*q.ind_welfare[1:q.aVP ,1:q.zP_low.-1])/sum(sim_0.measure[1:q.aVP ,1:q.zP_low.-1]),

welfare_avg_medR =  sum(sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1].*q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1])/sum(sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1]),
welfare_avg_medP =  sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].*q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1])/sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1]),
welfare_avg_medVP =  sum(sim_0.measure[1:q.aVP ,q.zP_low:q.zP_med.-1].*q.ind_welfare[1:q.aVP ,q.zP_low:q.zP_med.-1])/sum(sim_0.measure[1:q.aVP ,q.zP_low:q.zP_med.-1]),

welfare_avg_highR =  sum(sim_0.measure[q.aP:end,q.zP_med:end].*q.ind_welfare[q.aP:end,q.zP_med:end])/sum(sim_0.measure[q.aP:end,q.zP_med:end]),
welfare_avg_highP =  sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end].*q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end])/sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end]),
welfare_avg_highVP =  sum(sim_0.measure[1:q.aVP ,q.zP_med:end].*q.ind_welfare[1:q.aVP ,q.zP_med:end])/sum(sim_0.measure[1:q.aVP ,q.zP_med:end])))

# share of losers
q=merge(q,(share_losers_lowR = sum((q.ind_welfare[q.aP:end,1:q.zP_low.-1].<0) .* sim_0.measure[q.aP:end,1:q.zP_low.-1]),
share_losers_lowP = sum((q.ind_welfare[q.aVP.+1:q.aP.-1,1:q.zP_low.-1].<0) .* sim_0.measure[q.aVP.+1:q.aP.-1,1:q.zP_low.-1]),
share_losers_lowVP = sum((q.ind_welfare[1:q.aVP,1:q.zP_low.-1].<0) .*  sim_0.measure[1:q.aVP,1:q.zP_low-.1]),

share_losers_medR = sum((q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1].<0) .*  sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1]),
share_losers_medP = sum((q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].<0) .*sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1]),
share_losers_medVP = sum((q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1].<0) .*sim_0.measure[1:q.aVP,q.zP_low:q.zP_med.-1]),

share_losers_highR = sum((q.ind_welfare[q.aP:end,q.zP_med:end].<0) .*  sim_0.measure[q.aP:end,q.zP_med:end]),
share_losers_highP = sum((q.ind_welfare[q.aVP+1:q.aP-1,q.zP_med:end].<0) .* sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end]),
share_losers_highVP = sum((q.ind_welfare[1:q.aVP,q.zP_med:end].<0) .* sim_0.measure[1:q.aVP,q.zP_med:end])))

# share of winners
q=merge(q,(share_winners_lowR = sum((q.ind_welfare[q.aP:end,1:q.zP_low.-1].>=0) .* sim_0.measure[q.aP:end,1:q.zP_low.-1]),
share_winners_lowP = sum((q.ind_welfare[q.aVP.+1:q.aP.-1,1:q.zP_low.-1].>=0) .* sim_0.measure[q.aVP.+1:q.aP.-1,1:q.zP_low.-1]),
share_winners_lowVP = sum((q.ind_welfare[1:q.aVP,1:q.zP_low.-1].>=0) .*  sim_0.measure[1:q.aVP,1:q.zP_low.-1]),

share_winners_medR = sum((q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1].>=0) .*  sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1]),
share_winners_medP = sum((q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].>=0) .*sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1]),
share_winners_medVP = sum((q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1].>=0) .*sim_0.measure[1:q.aVP,q.zP_low:q.zP_med.-1]),

share_winners_highR = sum((q.ind_welfare[q.aP:end,q.zP_med:end].>=0) .*  sim_0.measure[q.aP:end,q.zP_med:end]),
share_winners_highP = sum((q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end].>=0) .* sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end]),
share_winners_highVP = sum((q.ind_welfare[1:q.aVP,q.zP_med:end].>=0) .* sim_0.measure[1:q.aVP,q.zP_med:end])))

# avg welfare of losers
q=merge(q,(welfare_avg_losers_lowR = sum(sim_0.measure[q.aP:end,1:q.zP_low-1] .* q.ind_welfare[q.aP:end,1:q.zP_low-1] .* (q.ind_welfare[q.aP:end,1:q.zP_low-1].<0))/sum(sim_0.measure[q.aP:end,1:q.zP_low-1].*(q.ind_welfare[q.aP:end,1:q.zP_low-1].<0)),
welfare_avg_losers_lowP =sum(sim_0.measure[q.aVP+1:q.aP-1,1:q.zP_low-1] .* q.ind_welfare[q.aVP+1:q.aP-1,1:q.zP_low-1] .* (q.ind_welfare[q.aVP+1:q.aP-1,1:q.zP_low-1].<0))/sum(sim_0.measure[q.aVP+1:q.aP-1,1:q.zP_low-1].*(q.ind_welfare[q.aVP+1:q.aP-1,1:q.zP_low-1].<0)),
welfare_avg_losers_lowVP = sum(sim_0.measure[1:q.aVP,1:q.zP_low-1] .* q.ind_welfare[1:q.aVP,1:q.zP_low-1] .* (q.ind_welfare[1:q.aVP,1:q.zP_low-1].<0))/sum(sim_0.measure[1:q.aVP,1:q.zP_low-1].*(q.ind_welfare[1:q.aVP,1:q.zP_low-1].<0)),

welfare_avg_losers_medR = sum(sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1] .* q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1] .* (q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1].<0))/sum(sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1].*(q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1].<0)),
welfare_avg_losers_medP = sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1] .* q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1] .* (q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].<0))/sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].*(q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].<0)),
welfare_avg_losers_medVP = sum(sim_0.measure[1:q.aVP,q.zP_low:q.zP_med.-1] .* q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1] .* (q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1].<0))/sum(sim_0.measure[1:q.aVP,q.zP_low:q.zP_med.-1].*(q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1].<0)),

welfare_avg_losers_highR = sum(sim_0.measure[q.aP:end,q.zP_med:end] .* q.ind_welfare[q.aP:end,q.zP_med:end] .* (q.ind_welfare[q.aP:end,q.zP_med:end].<0))/sum(sim_0.measure[q.aP:end,q.zP_med:end].*(q.ind_welfare[q.aP:end,q.zP_med:end].<0)),
welfare_avg_losers_highP = sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end] .* q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end] .* (q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end].<0))/sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end].*(q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end].<0)),
welfare_avg_losers_highVP = sum(sim_0.measure[1:q.aVP,q.zP_med:end] .* q.ind_welfare[1:q.aVP,q.zP_med:end] .* (q.ind_welfare[1:q.aVP,q.zP_med:end].<0))/sum(sim_0.measure[1:q.aVP,q.zP_med:end].*(q.ind_welfare[1:q.aVP,q.zP_med:end].<0))))

# avg welfare of winners
q=merge(q,(welfare_avg_winners_lowR = sum(sim_0.measure[q.aP:end,1:q.zP_low.-1] .* q.ind_welfare[q.aP:end,1:q.zP_low.-1] .* (q.ind_welfare[q.aP:end,1:q.zP_low.-1].>=0))/sum(sim_0.measure[q.aP:end,1:q.zP_low.-1].*(q.ind_welfare[q.aP:end,1:q.zP_low.-1].>=0)),
welfare_avg_winners_lowP = sum(sim_0.measure[q.aVP.+1:q.aP.-1,1:q.zP_low.-1] .* q.ind_welfare[q.aVP.+1:q.aP.-1,1:q.zP_low.-1] .* (q.ind_welfare[q.aVP.+1:q.aP.-1,1:q.zP_low.-1].>=0))/sum(sim_0.measure[q.aVP.+1:q.aP.-1,1:q.zP_low.-1].*(q.ind_welfare[q.aVP.+1:q.aP.-1,1:q.zP_low.-1].>=0)),
welfare_avg_winners_lowVP = sum(sim_0.measure[1:q.aVP,1:q.zP_low.-1] .* q.ind_welfare[1:q.aVP,1:q.zP_low.-1] .* (q.ind_welfare[1:q.aVP,1:q.zP_low.-1].>=0))/sum(sim_0.measure[1:q.aVP,1:q.zP_low.-1].*(q.ind_welfare[1:q.aVP,1:q.zP_low.-1].>=0)),

welfare_avg_winners_medR = sum(sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1] .* q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1] .* (q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1].>=0))/sum(sim_0.measure[q.aP:end,q.zP_low:q.zP_med.-1].*(q.ind_welfare[q.aP:end,q.zP_low:q.zP_med.-1].>=0)),
welfare_avg_winners_medP = sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1] .* q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1] .* (q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].>=0))/sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].*(q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_low:q.zP_med.-1].>=0)),
welfare_avg_winners_medVP = sum(sim_0.measure[1:q.aVP,q.zP_low:q.zP_med.-1] .* q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1] .* (q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1].>=0))/sum(sim_0.measure[1:q.aVP,q.zP_low:q.zP_med.-1].*(q.ind_welfare[1:q.aVP,q.zP_low:q.zP_med.-1].>=0)),

welfare_avg_winners_highR = sum(sim_0.measure[q.aP:end,q.zP_med:end] .* q.ind_welfare[q.aP:end,q.zP_med:end] .* (q.ind_welfare[q.aP:end,q.zP_med:end].>=0))/sum(sim_0.measure[q.aP:end,q.zP_med:end].*(q.ind_welfare[q.aP:end,q.zP_med:end].>=0)),
welfare_avg_winners_highP = sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end] .* q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end] .* (q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end].>=0))/sum(sim_0.measure[q.aVP.+1:q.aP.-1,q.zP_med:end].*(q.ind_welfare[q.aVP.+1:q.aP.-1,q.zP_med:end].>=0)),
welfare_avg_winners_highVP = sum(sim_0.measure[1:q.aVP,q.zP_med:end] .* q.ind_welfare[1:q.aVP,q.zP_med:end] .* (q.ind_welfare[1:q.aVP,q.zP_med:end].>=0))/sum(sim_0.measure[1:q.aVP,q.zP_med:end].*(q.ind_welfare[1:q.aVP,q.zP_med:end].>=0))))

q=merge(q,(table_exp=zeros(6,6),table_prod=zeros(6,6)))

### TABLES EXPORT STATUS & NET WORTH
q.table_exp[1,:] = [q.share_nonexpVP q.share_nonexpP q.share_nonexpR q.share_expVP q.share_expP q.share_expR]
q.table_exp[2,:] = [q.welfare_avg_nonexpVP q.welfare_avg_nonexpP q.welfare_avg_nonexpR q.welfare_avg_expVP q.welfare_avg_expP q.welfare_avg_expR]
q.table_exp[3,:] = [q.share_winners_nonexpVP q.share_winners_nonexpP q.share_winners_nonexpR q.share_winners_expVP q.share_winners_expP q.share_winners_expR]
q.table_exp[4,:] = [q.welfare_avg_winners_nonexpVP q.welfare_avg_winners_nonexpP q.welfare_avg_winners_nonexpR q.welfare_avg_winners_expVP q.welfare_avg_winners_expP q.welfare_avg_winners_expR]
q.table_exp[5,:] = [q.share_losers_nonexpVP q.share_losers_nonexpP q.share_losers_nonexpR q.share_losers_expVP q.share_losers_expP q.share_losers_expR]
q.table_exp[6,:] = [q.welfare_avg_losers_nonexpVP q.welfare_avg_losers_nonexpP q.welfare_avg_losers_nonexpR q.welfare_avg_losers_expVP q.welfare_avg_losers_expP q.welfare_avg_losers_expR]

### TABLES PRODUCTIVITY & NET WORTH

q.table_prod[1,:] = [q.share_lowVP q.share_lowP q.share_lowR q.share_medVP q.share_medP q.share_medR q.share_highVP q.share_highP q.share_highR]
q.table_prod[2,:] = [q.welfare_avg_lowVP q.welfare_avg_lowP q.welfare_avg_lowR q.welfare_avg_medVP q.welfare_avg_medP q.welfare_avg_medR q.welfare_avg_highVP q.welfare_avg_highP q.welfare_avg_highR]
q.table_prod[3,:] = [q.share_winners_lowVP q.share_winners_lowP q.share_winners_lowR q.share_winners_medVP q.share_winners_medP q.share_winners_medR q.share_winners_highVP q.share_winners_highP q.share_winners_highR]
q.table_prod[4,:] = [q.welfare_avg_winners_lowVP q.welfare_avg_winners_lowP q.welfare_avg_winners_lowR q.welfare_avg_winners_medVP q.welfare_avg_winners_medP q.welfare_avg_winners_medR q.welfare_avg_winners_highVP q.welfare_avg_winners_highP q.welfare_avg_winners_highR]
q.table_prod[5,:] = [q.share_losers_lowVP q.share_losers_lowP q.share_losers_lowR q.share_losers_medVP q.share_losers_medP q.share_losers_medR q.share_losers_highVP q.share_losers_highP q.share_losers_highR]
q.table_prod[6,:] = [q.welfare_avg_losers_lowVP q.welfare_avg_losers_lowP q.welfare_avg_losers_lowR q.welfare_avg_losers_medVP q.welfare_avg_losers_medP q.welfare_avg_losers_medR q.welfare_avg_losers_highVP q.welfare_avg_losers_highP q.welfare_avg_losers_highR]

### Bins
sav = 0
zbin = 5
abin = 5

q=merge(q,(measure_bins_5=zeros(abin,zbin),measure_bins_5_end=zeros(abin,zbin),avg_welfare_bins_5=zeros(abin,zbin),winners_bins_5=zeros(abin,zbin),losers_bins_5=zeros(abin,zbin)))

for j = 1:zbin
    for i = 1:abin
        pL = (j-1)/zbin*s.z_grid_size+1
        pH = j/zbin*s.z_grid_size
        aL = (i-1)/abin*s.a_grid_size+1
        aH = i/abin*s.a_grid_size

        q.measure_bins_5[i,j] = sum(sim_0.measure[aL:aH,pL:pH]);
        q.measure_bins_5_end[i,j] = sum(sim_end.measure[aL:aH,pL:pH]),
        q.avg_welfare_bins_5[i,j] = sum(sim_0.measure[aL:aH,pL:pH] .* q.ind_welfare[aL:aH,pL:pH])./sum(sim_0.measure[aL:aH,pL:pH]),
        q.winners_bins_5[i,j] = sum(sim_0.measure[aL:aH,pL:pH] .* (q.ind_welfare[aL:aH,pL:pH].>=0))/sum(sim_0.measure[aL:aH,pL:pH]),
        q.losers_bins_5[i,j] = sum(sim_0.measure[aL:aH,pL:pH] .* (q.ind_welfare[aL:aH,pL:pH].<0))/sum(sim_0.measure[aL:aH,pL:pH])
    end
end

#FROM HERE UNTIL "3D bar (10 bins)" PENDING

figure(1)
b1 = bar3(q.measure_bins_5); colorbar
for k = 1:length(b1)
    zdata = b1(k).ZData;
    b1(k).CData = zdata;
    b1(k).FaceColor = 'interp';
end
ylabel('assets')
xlabel('productivity')
# view(32, 22)
if sav == 1
        saveas(gcf,'Figures/measure_FF_tau_c_5','epsc')
        saveas(gcf,'Figures/measure_FF_tau_c_5','pdf')
end

figure(2)
b2 = bar3(q.avg_welfare_bins_5); colorbar
for k = 1:length(b2)
    zdata = b2(k).ZData;
    b2(k).CData = zdata;
    b2(k).FaceColor = 'interp';
end
xlabel('productivity')
ylabel('assets')
# view(32, 22)
if sav == 1
        saveas(gcf,'Figures/welfare_FF_tau_c_5','epsc')
        saveas(gcf,'Figures/welfare_FF_tau_c_5','pdf')
end

figure(3)
b3 = bar3(q.losers_bins_5); colorbar
for k = 1:length(b3)
    zdata = b3(k).ZData;
    b3(k).CData = zdata;
    b3(k).FaceColor = 'interp';
end
xlabel('productivity')
ylabel('assets')
set(gca,'Zlim',[0,1])
caxis([0 1])
# view(32, 22)
if sav == 1
        saveas(gcf,'Figures/welfare_FF_tau_c_5','epsc')
        saveas(gcf,'Figures/welfare_FF_tau_c_5','pdf')
end

figure(4)
imagesc(q.winners_bins_5); colorbar
xlabel('productivity')
ylabel('assets')
# view(32, 22)
if sav == 1
        saveas(gcf,'Figures/winners_FF_tau_c_5','epsc')
        saveas(gcf,'Figures/winners_FF_tau_c_5','pdf')
end

figure(5)
imagesc(q.ind_welfare); colorbar
if sav == 1
        saveas(gcf,'Figures/welfare_image','epsc')
        saveas(gcf,'Figures/welfare_image','pdf')
end

### 3D bar (10 bins)
zbin = 10;
abin = 10;

q=merge(q,(measure_bins_10=zeros(abin,zbin),avg_welfare_bins_10=zeros(abin,zbin),winners_bins_10=zeros(abin,zbin),losers_bins_10=zeros(abin,zbin)))

for j = 1:zbin
    for i = 1:abin
        pL = (j-1)/zbin*s.z_grid_size+1;
        pH = j/zbin*s.z_grid_size;
        aL = (i-1)/abin*s.a_grid_size+1;
        aH = i/abin*s.a_grid_size;

        q.measure_bins_10[i,j] = sum(sim_0.measure[aL:aH,pL:pH])
        q.avg_welfare_bins_10[i,j] = sum(sim_0.measure[aL:aH,pL:pH] .* q.ind_welfare[aL:aH,pL:pH])./sum(sim_0.measure[aL:aH,pL:pH])
        q.winners_bins_10[i,j] = sum(sim_0.measure[aL:aH,pL:pH] .* (q.ind_welfare[aL:aH,pL:pH].>=0))/ sum(sim_0.measure[aL:aH,pL:pH])
        q.losers_bins_10[i,j] = sum(sim_0.measure[aL:aH,pL:pH] .* (q.ind_welfare[aL:aH,pL:pH].<0))/sum(sim_0.measure[aL:aH,pL:pH])
    end
end
if s.graphs_matlab ==1
file=matopen("data_welfaregraphs.mat","w")
write(file,"measure_bins_10",q.measure_bins_10)
close(file)
else

#FROM HERE UNTIL THE END PENDING
figure(5)
b1 = bar3(q.measure_bins_10); colorbar
for k = 1:length(b1)
    zdata = b1(k).ZData;
    b1(k).CData = zdata;
    b1(k).FaceColor = 'interp';
end
ylabel('assets')
xlabel('productivity')
# view(32, 22)
if sav == 1
        saveas(gcf,'Figures/measure_FF_tau_c_10','epsc')
        saveas(gcf,'Figures/measure_FF_tau_c_10','pdf')
end

figure(6)
b2 = bar3(q.avg_welfare_bins_10); colorbar
for k = 1:length(b2)
    zdata = b2(k).ZData;
    b2(k).CData = zdata;
    b2(k).FaceColor = 'interp';
end
xlabel('productivity')
ylabel('assets')
# view(32, 22)
if sav == 1
        saveas(gcf,'Figures/welfare_FF_tau_c_10','epsc')
        saveas(gcf,'Figures/welfare_FF_tau_c_10','pdf')
end

figure(7)
imagesc(q.winners_bins_10); colorbar
# xlabel('productivity')
# ylabel('assets')
# view(32, 22)
if sav == 1
        saveas(gcf,'Figures/winners_FF_tau_c_10','epsc')
        saveas(gcf,'Figures/winners_FF_tau_c_10','pdf')
end

end


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

JLLWrappers.@generate_main_file_header("Cbc")
JLLWrappers.@generate_main_file("Cbc", UUID("38041ee0-ae04-5750-a4d2-bb4d0d83d27d"))
end  # module Cbc_jll


text:
### Compute empirical covariance matrix recursively

covariance(lastC::Real, k::Integer, x::Real, lastmean::Real, secondlastmean::Real) =
  ((k-1)*lastC+abs2(x)-(k+1)*abs2(lastmean)+k*abs2(secondlastmean))/k

function covariance!(
  C::RealMatrix,
  lastC::RealMatrix,
  k::Integer,
  x::RealVector,
  lastmean::RealVector,
  secondlastmean::RealVector
)
  C[:, :] = (k-1)*lastC
  BLAS.ger!(1.0, x, x, C)
  BLAS.ger!(1.0, -(k+1)*lastmean, lastmean, C)
  BLAS.ger!(1.0, k*secondlastmean, secondlastmean, C)
  C[:, :] /= k
end


text:
using LPVSpectral, Statistics, Random
export Sinusoidal

"""
    Sinusoidal(fs)
Decompose a timeseries `s` into a **sum** `x + r`, where `x`
are sinusoidal components with the given frequencies `fs` that minimize
coefficients ``A, \\phi`` of the expression
```math
s \\approx A_0 + \\sum_i A_i \\cos(2\\pi f_i t + \\phi_i)
```
with ``A_0 = \\bar{s}`` the mean.

This method uses a new least-squares algorithm in frequency domain using the package
[LPVSpectral.jl](https://github.com/baggepinnen/LPVSpectral.jl), see[^Bagge2017].
It works for non-equispaced `t` axis (and also normal), is generally very accurate
(if choosen frequencies are not too close), but has performance scaling of
O(N^2.4) instead of O(n log(n)) of [`Fourier`](@ref).

Because it can work with arbitrary signal length the method always estimates the
zero-frequency Fourier component, and attributes it to `x`.
The fitted coefficients ``A, \\phi`` are available as fields `.A` and `.φ` of the struct
(first entry is zero-frequency component ``A_0, \\phi_0``).

[^Bagge2017]: F. Bagge Carlson et al., [Linear Parameter-Varying Spectral Decomposition](https://lup.lub.lu.se/search/publication/ac32368e-e199-44ff-b76a-36668ac7d595).
"""
struct Sinusoidal{F<:AbstractArray{<:AbstractFloat}} <: Decomposition
    fs::F
    A::Vector{Float64}
    φ::Vector{Float64}
end
Sinusoidal(fs) = Sinusoidal(fs, zeros(length(fs)+1), zeros(length(fs)+1))

function decompose(t, s, method::Sinusoidal)
    @assert 0 ∉ method.fs
    fs = [0, method.fs...]
    χ, = ls_spectral(s, t, fs)
    As = abs.(χ) ./ sqrt(2*length(fs))
    φs = angle.(χ)
    method.A .= As
    method.φ .= φs
    x = fill(As[1], length(s))
    for i in 2:length(fs)
        @. x += As[i]*cos(2π*fs[i]*t + φs[i])
    end
    return x, s .- x
end


text:
println("********** FWI::testSensitivity ****************")
import FactoredEikonalFastMarching.getAnalyticalConstGradInv2D
function checkDerivativeFWI(f::Function,x0,param,Iact;out::Bool=true,tol::Number=1.9,nSuccess::Int=3)
	if out
		println(@sprintf("%9s\t%9s\t%9s\t%9s\t%9s\t%5s","h","E0","E1","O1","O2","OK?"))
	end
	v = 5.0*randn(size(x0))
	v = Iact*Iact'*v;
	# if eltype(x0)<:Complex
		# v += 1im*randn(size(v))
	# end
	f0,dvf = f(x0,v,param)

	# dvf   = real(dvf)
	Error = zeros(10,2)
	Order = zeros(10,2)
	Success = zeros(10)
	for j=1:6
		ft = f(x0.+(10.0^(-j))*v,zeros(size(x0)),param)[1]            # function value
		Error[j,1] = norm(f0.-ft)/norm(f0)           # Error TaylorPoly 0
		Error[j,2] = norm(f0 .+(10.0^(-j))*dvf .- ft)/norm(f0) # Error TaylorPoly 1
		if j>1
			Order[j,:] .= log10.(Error[j-1,:]./Error[j,:]);
		end
		if (Order[j,2]>tol) || (Error[j,1]/Error[j,2] > 100); Success[j]=1; end
		if out
			println(@sprintf("%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%5d",
							10.0^(-j), Error[j,1],Error[j,2], Order[j,1],Order[j,2],Success[j]))
		end
	end
	pass = sum(Success) > nSuccess
	return  pass,Error,Order
end

println("########### 2D TEST SETTING ##############################");
Omega = [0.0,4.0,0.0,8.0];
Minv = getRegularMesh(Omega,[64,128]);
println(Minv.n)
m = getAnalyticalConstGradInv2D(Minv.n,Minv.h)[1];
m = m';

pad = 16;
jump = 15;
offset = 100;

# omega = [0.5./maximum(Minv.h), 0.2./maximum(Minv.h)];
omega = [0.5./maximum(Minv.h)];

WaveletCoef = rand(ComplexF64,size(omega));
println(omega)

RCVfile = string("rcvMap.dat");
SRCfile = string("srcMap.dat");
writeSrcRcvLocFile(SRCfile,Minv,pad,jump);
writeSrcRcvLocFile(RCVfile,Minv,pad,2);
srcNodeMap = readSrcRcvLocationFile(SRCfile,Minv);
rcvNodeMap = readSrcRcvLocationFile(RCVfile,Minv);
Q = generateSrcRcvProjOperators(Minv.n.+1,srcNodeMap);
P = generateSrcRcvProjOperators(Minv.n.+1,rcvNodeMap);

ABLPad = pad.+5;

Minv.n.-=1;
gamma = Helmholtz.getABL(Minv,true,ones(Int64,Minv.dim)*ABLPad,omega[1]);
Minv.n.+=1;

attenuation = 0.01*omega[1];
gamma .+= attenuation; # adding Attenuation.

# Ainv = getMUMPSsolver([],0,0,0);
Ainv = getJuliaSolver();

###################################################################################################
# levels      = 2;
# numCores 	= 4;
# BLAS.set_num_threads(numCores);
# maxIter     = 50;
# relativeTol = 1e-5;
# relaxType   = "SPAI";
# relaxParam  = 1.0;
# relaxPre 	= 2;
# relaxPost   = 2;
# cycleType   ='W';
# coarseSolveType = "NoMUMPS";
# MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType,0.0,0.0);
# shift = 0.15;
# Hparam = HelmholtzParam(Minv,zeros(0),zeros(0),0.0,true,true);
# Ainv = getShiftedLaplacianMultigridSolver(Hparam, MG,shift);
###################################################################################################

println("Number of Sources: ",size(Q,2));

(pFor,contDiv,SourcesSubInd) = getFWIparam(omega,WaveletCoef,vec(gamma),Q,P,Minv,Ainv,workers(),16);

# setup active cells - this is important so that Sommerfeld BC will not interfere with the sensitivities.

N = prod(Minv.n);
Iact = sparse(1.0I, N, N);
mask = zeros(N);
mask[gamma[:] .>= 0.2*maximum(gamma)] .= 1.0;
Iact = Iact[:,mask .== 0.0];

function f(m,v=zeros(size(m)),pFor = [])
	(D,pFor) = getData(m[:],pFor);
	D = fetch(D[1]);
	pf = fetch(pFor[1]);
	Jv = 0;
	if norm(v)>0.0
		Jv = getSensMatVec(v[:],m[:],pf);
	end
	return D[:],Jv;
end
checkDerivativeFWI(f,m[:],pFor,Iact);

pf = fetch(pFor[1]);
v = rand(N);
Jv = getSensMatVec(v,m[:],pf);
Z = rand(ComplexF64,size(Jv));
JtZ = getSensTMatVec(Z,m[:],pf);

println(dot(Jv[:],Z[:]))
println(dot(JtZ[:],v))

@test real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) < 1e-5
if real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) >= 1e-5
	println("Checking that (J^T)^T = J for FWI falied: ",real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) );
end

############################## STOCHASTIC #####################################################
# println("####################### STOCHASTIC SELECTION ################################")
# sourcesSelection = randperm(size(Q,2))[1:div(size(Q,2),2)];
# setSourceSelection(pFor[1], sourcesSelection)

# function f(m,v=zeros(size(m)),pFor = [])
	# (D,pFor) = getData(m[:],pFor,Mesh2Mesh);
	# D = fetch(D[1]);
	# pf = fetch(pFor[1]);
	# Jv = 0;
	# if norm(v)>0.0
		# Jv = getSensMatVec(v[:],m[:],pf);
	# end
	# return D[:],Jv;
# end
# checkDerivativeFWI(f,m[:],pFor,Iact);

# pf = fetch(pFor[1]);
# v = rand(N);
# Jv = getSensMatVec(v,m[:],pf);
# Z = rand(Complex128,size(Jv));
# JtZ = getSensTMatVec(Z,m[:],pf);

# println(dot(Jv[:],Z[:]))
# println(dot(JtZ[:],v))

# if real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) > 1e-8
	# warn("Checking that (J^T)^T = J for FWI falied: ",real(dot(JtZ[:],v) - dot(Jv[:],Z[:])) );
# end

rm("srcMap.dat");
rm("rcvMap.dat");


text:
export GoalReward, BonusGoalReward

"Specification where reaching the goal delivers reward."
@kwdef struct GoalReward <: Goal
    terms::Vector{Term} # Goal terms to be satisfied
    reward::Float64 = 1.0 # Reward gained from reaching goal
    discount::Float64 = 0.9 # Discount factor
end

GoalReward(problem::Problem) = GoalReward(flatten_conjs(problem.goal), 1.0, 0.9)
GoalReward(terms, reward) = GoalReward(terms, reward, 0.9)
GoalReward(terms) = GoalReward(terms, 1.0, 0.9)
GoalReward(term::Term, reward, discount) =
    GoalReward(flatten_conjs(term), reward, discount)

Base.hash(spec::GoalReward, h::UInt) =
    hash(spec.reward, hash(spec.discount, hash(Set(spec.terms), h)))
Base.isequal(s1::GoalReward, s2::GoalReward) =
    s1.reward == s2.reward && s1.discount == s2.discount &&
    Set(s1.terms) == Set(s2.terms)

is_goal(spec::GoalReward, domain::Domain, state::State) =
    satisfy(spec.terms, state, domain)[1]
is_violated(spec::GoalReward, domain::Domain, state::State) = false
get_cost(spec::GoalReward, domain::Domain, s1::State, a::Term, s2::State) =
    -get_reward(spec, domain, s1, a, s2)
get_reward(spec::GoalReward, domain::Domain, s1::State, a::Term, s2::State) =
    is_goal(spec, domain, s2) ? spec.reward : 0.0
get_discount(spec::GoalReward) = spec.discount
get_goal_terms(spec::GoalReward) = spec.terms

discounted(spec::GoalReward, discount::Float64) =
    GoalReward(spec.terms, spec.reward, discount * spec.discount)

"Wrapper specification that delivers additional reward upon reaching goal."
@kwdef struct BonusGoalReward{G <: Goal} <: Goal
    goal::G # Goal specification to be satisfied
    reward::Float64 = 1.0 # Additional reward gained from reaching goal
    discount::Float64 = 1.0 # Discount factor
end

BonusGoalReward(spec::Goal) = BonusGoalReward(spec, 1.0, 1.0)
BonusGoalReward(spec::Goal, reward) = BonusGoalReward(spec, reward, 1.0)
BonusGoalReward(spec::BonusGoalReward, reward, discount) =
    BonusGoalReward(spec.goal, reward + spec.reward, discount * spec.discount)

is_goal(spec::BonusGoalReward, domain::Domain, state::State) =
    is_goal(spec.goal, domain, state)
is_violated(spec::BonusGoalReward, domain::Domain, state::State) =
    is_violated(spec.goal, domain, state)
get_cost(spec::BonusGoalReward, domain::Domain, s1::State, a::Term, s2::State) =
    -get_reward(spec, domain, s1, a, s2)
get_reward(spec::BonusGoalReward, domain::Domain, s1::State, a::Term, s2::State) =
    get_reward(spec.goal, domain, s1, a, s2) +
    is_goal(spec, domain, s2) ? spec.reward : 0.0
get_discount(spec::BonusGoalReward) = spec.discount
get_goal_terms(spec::BonusGoalReward) = get_goal_terms(spec.goal)

discounted(spec::BonusGoalReward, discount::Float64) =
    BonusGoalReward(spec.goal, spec.reward, discount * spec.discount)


text:
import Plots._current_plots_version

# replace `f(args...)` with `f(rng, args...)` for `f ∈ (rand, randn)`
function replace_rand!(ex) end
function replace_rand!(ex::Expr)
    for arg in ex.args
        replace_rand!(arg)
    end
    if ex.head === :call && ex.args[1] ∈ (:rand, :randn, :(Plots.fakedata))
        pushfirst!(ex.args, ex.args[1])
        ex.args[2] = :rng
    end
end
function fix_rand!(ex)
    replace_rand!(ex)
    pushfirst!(ex.args[1].args, :(rng = StableRNG(PLOTS_SEED)))
end

function image_comparison_tests(
    pkg::Symbol,
    idx::Int;
    debug = false,
    popup = !is_ci(),
    sigma = [1, 1],
    tol = 1e-2,
)
    Plots._debugMode.on = debug
    example = Plots._examples[idx]
    Plots.theme(:default)
    @info("Testing plot: $pkg:$idx:$(example.header)")
    backend(pkg)
    backend()
    default(size = (500, 300))

    fn = "ref$idx.png"
    reffn = reference_file(pkg, idx, _current_plots_version)
    newfn = joinpath(reference_path(pkg, _current_plots_version), fn)
    @debug example.exprs

    # test function
    func = (fn, idx) -> begin
        expr = Expr(:block)
        append!(expr.args, example.exprs)
        fix_rand!(expr)
        eval(expr)
        png(fn)
    end

    # the test
    vtest = VisualTest(func, reffn, idx)
    test_images(vtest, popup = popup, sigma = sigma, tol = tol, newfn = newfn)
end

function image_comparison_facts(
    pkg::Symbol;
    skip = [],      # skip these examples (int index)
    only = nothing, # limit to these examples (int index)
    debug = false,  # print debug information?
    sigma = [1, 1],  # number of pixels to "blur"
    tol = 1e-2,
)     # acceptable error (percent)
    for i in 1:length(Plots._examples)
        i in skip && continue
        if only === nothing || i in only
            @test image_comparison_tests(pkg, i, debug = debug, sigma = sigma, tol = tol) |>
                  success == true
        end
    end
end


text:
################
# Comment Type #
################

type Comment <: GitLabType
    noteable_type::Nullable{GitLabString}
    ## created_at::Nullable{Dates.DateTime}
    created_at::Nullable{GitLabString}
    line_code::Nullable{Int64}
    note::Nullable{GitLabString}
    author::Nullable{Owner}
    author_id::Nullable{Int64}
    updated_by_id::Nullable{GitLabString}
    noteable_id::Nullable{Int64}
    commit_id::Nullable{GitLabString}
    system::Nullable{Bool}
    url::Nullable{HttpCommon.URI}
    is_award::Nullable{Bool}
    st_diff::Nullable{GitLabString}
    id::Nullable{Int64}
    ## updated_at::Nullable{Dates.DateTime}
    updated_at::Nullable{GitLabString}
    project_id::Nullable{Int64}
    attachment::Nullable{GitLabString}
    ## type::Nullable{GitLabString}
    path::Nullable{GitLabString}
    line::Nullable{GitLabString}
    line_type::Nullable{GitLabString}

#=
    original_commit_id::Nullable{GitLabString}
    original_position::Nullable{Int}
    position::Nullable{Int}
    html_url::Nullable{HttpCommon.URI}
    issue_url::Nullable{HttpCommon.URI}
    pull_request_url::Nullable{HttpCommon.URI}
    user::Nullable{Owner}
=#

#=
    ## Issue comment
    {
        "note": "Test ...",
        "path": null,
        "line": null,
        "line_type": null,
        "author": {
          "name": "Pradeep",
          "username": "mdpradeep",
          "id": 2,
          "state": "active",
          "avatar_url": "http://www.gravatar.com/avatar/7e32a35a20817e0258e12665c9099422?s=80&d=identicon",
          "web_url": "http://104.197.141.88/u/mdpradeep"
        },
        "created_at": "2016-07-16T16:02:12.923Z"
    }

    ## Merge request comment
  "noteable_type" => "MergeRequest"
  "created_at"    => "2016-07-19 09:45:27 UTC"
  "line_code"     => nothing
  "note"          => "`sayhello(\"ABC\", \"cool\")`"
  "author_id"     => 2
  "updated_by_id" => nothing
  "noteable_id"   => 4
  "commit_id"     => ""
  "system"        => false
  "url"           => "http://104.197.141.88/mdpradeep/TestProject1/merge_requests/4#note_134"
  "is_award"      => false
  "st_diff"       => nothing
  "id"            => 134
  "updated_at"    => "2016-07-19 09:45:27 UTC"
  "project_id"    => 1
  "attachment"    => nothing
  "type"          => nothing

    ## commit comment
  "noteable_type" => "Commit"
  "created_at"    => "2016-07-18 06:49:53 UTC"
  "line_code"     => nothing
  "note"          => "`sayhello(\\\"ABC\\\", \\\"cool\\\")`"
  "author_id"     => 2
  "updated_by_id" => nothing
  "noteable_id"   => nothing
  "commit_id"     => "d1d585d2fbc0c3a052e219f30504be8d1621a2a9"
  "system"        => false
  "url"           => "http://104.197.141.88/mdpradeep/TestProject1/commit/d1d585d2fbc0c3a052e219f30504be8d1621a2a9#note_46"
  "is_award"      => false
  "st_diff"       => nothing
  "id"            => 46
  "updated_at"    => "2016-07-18 06:49:53 UTC"
  "project_id"    => 1
  "attachment"    => nothing
  "type"          => nothing
=#

end

Comment(data::Dict) = json2gitlab(Comment, data)
Comment(id::Int) = Comment(Dict("id" => id))

namefield(comment::Comment) = comment.id

kind_err_str(kind) = ("Error building comment request: :$kind is not a valid kind of comment.\n"*
                      "The only valid comment kinds are: :issue, :review, :commit")
###############
# API Methods #
###############

function comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        ## MDP path = "/repos/$(name(repo))/issues/comments/$(name(item))"
        path = "/api/v3/projects/$(get(repo.id))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/api/v3/projects/$(get(repo.id))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/api/v3/projects/$(get(repo.id))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_get_json(path; options...))
end

function comments(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        ## MDP path = "/repos/$(name(repo))/issues/$(name(item))/comments"
        path = "/api/v3/projects/$(get(repo.id))/issues/$(name(item))/notes"
    elseif kind == :review
        ## MDP path = "/api/v3/projects/$(get(repo.id))/pulls/$(name(item))/comments"
        path = "/api/v3/projects/$(get(repo.id))/merge_requests/$(name(item))/notes"
    elseif kind == :commit
        path = "/api/v3/projects/$(get(repo.id))/repository/commits/$(name(item))/comments"
    else
        error(kind_err_str(kind))
    end

    results, page_data = gh_get_paged_json(path; options...)
    return map(Comment, results), page_data
end

function create_comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        ## MDP path = "/repos/$(name(repo))/issues/$(name(item))/comments"
        path = "/api/v3/projects/$(get(repo.id))/issues/$(name(item))/notes"
    elseif kind == :review
        ## MDP path = "/repos/$(name(repo))/pulls/$(name(item))/comments"
        path = "/api/v3/projects/$(get(repo.id))/merge_requests/$(name(item))/notes"
    elseif kind == :commit
        ## MDP path = "/repos/$(name(repo))/commits/$(name(item))/comments"
        path = "/api/v3/projects/$(get(repo.id))/repository/commits/$(name(item))/comments"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_post_json(path; options...))
end

function edit_comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        ## MDP path = "/repos/$(name(repo))/issues/comments/$(name(item))"
        path = "/api/v3/projects/$(get(repo.id))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/api/v3/projects/$(get(repo.id))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/api/v3/projects/$(get(repo.id))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return Comment(gh_patch_json(path; options...))
end

function delete_comment(repo, item, kind = :issue; options...)
    if (kind == :issue) || (kind == :pr)
        path = "/api/v3/projects/$(get(repo.id))/issues/comments/$(name(item))"
    elseif kind == :review
        path = "/api/v3/projects/$(get(repo.id))/pulls/comments/$(name(item))"
    elseif kind == :commit
        path = "/api/v3/projects/$(get(repo.id))/comments/$(name(item))"
    else
        error(kind_err_str(kind))
    end
    return gh_delete(path; options...)
end


text:
const VFTOL = 1e-10

"""
    vfit!(EV0, t, ddm)

Given `ubV(t) ≡ u + βV(x')`, update `EV0 ← E[max u + β V(x')]`
"""
function vfit!(EV0, t::DCDPTmpVars, ddm::DynamicDrillModel)
    if anticipate_t1ev(ddm)
        logsumexp_and_softmax!(lse(t), q(t), tmp(t), ubV(t))
    else
        maximum!(add_1_dim(lse(t)), ubV(t))
    end
    A_mul_B_md!(EV0, ztransition(ddm), lse(t), 1)
end

# preserves ubV & updates derivatives
function vfit!(EV0, dEV0, t::DCDPTmpVars, ddm::DynamicDrillModel)
    if anticipate_t1ev(ddm)
        logsumexp_and_softmax!(lse(t), q(t), tmp(t), ubV(t))
    else
        findmax!( add_1_dim(lse(t)), add_1_dim(tmp_cart(t)), ubV(t) )
        fill!( q(t), 0)
        @inbounds @simd for i in tmp_cart(t)
            setindex!(q(t), 1, i)
        end
    end
    A_mul_B_md!(EV0, ztransition(ddm), lse(t), 1)

    sumdubV = view(dubVperm(t), :,:,:,1)
    sumprod!(sumdubV, dubVperm(t), q(t))
    A_mul_B_md!(dEV0, ztransition(ddm), sumdubV, 1)
end

# --------------------------- VFIT until conv ----------------------------

function solve_inf_vfit!(EV0, t::DCDPTmpVars, ddm::DynamicDrillModel; maxit=60, vftol=VFTOL, kwargs...)

    iter = zero(maxit)
    while true
        vfit!(tmp(t), t, ddm)
        bnds = extrema( tmp(t) .- EV0 ) .* beta_1minusbeta(ddm)
        ubV(t)[:,:,1] .= discount(ddm) .* (EV0 .= tmp(t))
        iter += 1
        converged = all(abs.(bnds) .< vftol)
        if converged  ||  iter >= maxit
            return converged, iter, bnds
        end
    end

end

# --------------------------- basic pfit ----------------------------

function div_and_update_direct!(EV0, t, ddm)
    J = size(EV0, 2)
    ΔEV = lse(t)
    q0 = ubV(t)

    for j in OneTo(J)
        q0j  = uview(q0, :, j, 1)
        ΔEVj = uview(ΔEV, :, j)

        update_IminusTVp!(t, ddm, q0j)
        fact = lu(IminusTEVp(t))
        ldiv!(fact, ΔEVj)  # Vtmp = [I - T'(V)] \ [V - T(V)]
    end
    EV0 .-= ΔEV # update V
end

function div_and_update_indirect!(EV0, t, ddm, vftol; kwargs...)

    J = size(EV0,2)
    ΔEV = lse(t)
    q0 = ubV(t)
    ΔEVjnew = tmpEVj(t)

    for j in OneTo(J)
        q0j  = uview(q0, :, j, 1)
        ΔEVj = uview(ΔEV, :, j)
        if norm(ΔEVj, Inf) > vftol*1e-2  # b/c if ΔEVj ≈ 0, then we blow up
            fill!(ΔEVjnew, 0)
            update_IminusTVp!(t, ddm, q0j)
            # bicgstabl! results in NaNs
            gmres!(ΔEVjnew, IminusTEVp(t), ΔEVj; initially_zero=true, kwargs...)
            copyto!(ΔEVj, ΔEVjnew)
            # all(isfinite.(ΔEVjnew)) || throw(error("ΔEVjnew not finite. j=$j, ΔEVj = $ΔEVj, ΔEVjnew = $ΔEVjnew"))
        end
    end
    EV0 .-= ΔEV # update V
end

function pfit!(EV0::AbstractMatrix, t::DCDPTmpVars, ddm::DynamicDrillModel; vftol=VFTOL, kwargs...)

    ΔEV = lse(t)
    q0 = ubV(t)

    if anticipate_t1ev(ddm)
        logsumexp_and_softmax!(ΔEV, q0, tmp(t), ubV(t), 1)
    else
        findmax!(add_1_dim(ΔEV), add_1_dim(tmp_cart(t)), ubV(t))
        q0[:,:,1] .= last.(getfield.(tmp_cart(t), :I)) .== 1         # update q0 as Pr(d=0|x)
    end
    A_mul_B_md!(tmp(t), ztransition(ddm), ΔEV, 1)

    # compute difference & check bnds
    bnds = extrema(ΔEV .= EV0 .- tmp(t)) .* -beta_1minusbeta(ddm)
    if all(abs.(bnds) .< vftol)
        copyto!(EV0, tmp(t))
        return bnds
    end
    # Vtmp = [I - T'(V)] \ [V - T(V)]
    # V .= -Vtmp
    # div_and_update_direct!(EV0, t, ddm)
    div_and_update_indirect!(EV0, t, ddm, vftol; abstol=1e-13)
    return extrema(ΔEV) .* -beta_1minusbeta(ddm) # get norm
end

# --------------------------- pfit until convergence ----------------------------

function solve_inf_pfit!(EV0::AbstractMatrix, t::DCDPTmpVars, ddm::DynamicDrillModel; maxit=40, vftol=VFTOL, kwargs...)
    iter = zero(maxit)
    while true
        bnds = pfit!(EV0, t, ddm; vftol=vftol)
        iter += 1
        converged = all(abs.(bnds) .< vftol)
        if converged  ||  iter >= maxit
            return converged, iter, bnds
        end
        ubVfull(t)[:,:,1] .= discount(ddm) .* EV0
    end
end

# --------------------------- inf horizon gradient ----------------------------

# note -- this destroys ubV
function gradinf!(dEV0::AbstractArray3, t::DCDPTmpVars, ddm::DynamicDrillModel)

    nz, nψ, nk, nd = size(dubVperm(t))
    nd >= 2 || throw(error("Need dubV with at least 2+ action possibilities"))

    sumdubV   = view(dubVperm(t), :,:,:,1)
    ΠsumdubV  = view(dubVperm(t), :,:,:,2)

    if anticipate_t1ev(ddm)
        softmax3!(ubV(t), lse(t), tmp(t), ubV(t))
        # softmax3!(q, lse, tmp)
    else
        findmax!(add_1_dim(lse(t)), add_1_dim(tmp_cart(t)), ubV(t))
        fill!(ubV(t), 0)
        @inbounds for i in tmp_cart(t)
            setindex!(ubV(t), 1, i)
        end
    end

    # for dubV/dθt
    sumprod!(sumdubV, dubVperm(t), ubV(t))
    A_mul_B_md!(ΠsumdubV, ztransition(ddm), sumdubV, 1)
    gradinf_inner_direct!(dEV0, ΠsumdubV, t, ddm)
    # gradinf_inner_indirect!(dEV0, sumdubV, ΠsumdubV, t, ddm; gradinfmaxit=2)
end

function gradinf_inner_direct!(dEV0, ΠsumdubV, t, ddm)
    nz, nψ, nk = size(dEV0)

    for j in OneTo(nψ)
        qj = view(ubV(t), :, j, 1)
        update_IminusTVp!(t, ddm, qj)

        fact = lu(IminusTEVp(t))

        # Note: cannot do this with @view(dEV0[:,j,:])
        ΠsumdubVj(t) .= view(ΠsumdubV, :, j, :)
        ldiv!(dev0tmpj(t), fact, ΠsumdubVj(t)) # ΠsumdubV[:,j,:])
        dEV0[:,j,:] .= dev0tmpj(t)
    end
end

function gradinf_inner_indirect!(dEV0, sumdubV, ΠsumdubV, t, ddm; gradinfmaxit=1)
    nz, nψ, nθt = size(dEV0)

    # tmpvars
    βqdEV0 = sumdubV
    q0 = view(ubV(t), :, :, 1)

    # get some good starting values for iterative solver
    dEV0 .= ΠsumdubV
    for it = 1:gradinfmaxit
        βqdEV0 .= discount(ddm) .* q0 .* dEV0
        A_mul_B_md!(dEV0, ztransition(ddm), βqdEV0, 1)
        dEV0 .+= ΠsumdubV
    end

    # solve linear system for infinite horizon gradient
    for j in OneTo(nψ)
        qj = view(ubV(t), :, j, 1)
        update_IminusTVp!(t, ddm, qj)
        for k = OneTo(nθt)
            dEV0jk = uview(dEV0, :, j, k)
            ΠsumdubVjk = uview(ΠsumdubV, :, j, k)
            if any(dEV0jk != 0)
                gmres!(dEV0jk, IminusTEVp(t), ΠsumdubVjk; initially_zero=false)
            end
        end
    end
end

# # --------------------------- double vfit/pfit loop -----------------------

function solve_inf_vfit_pfit!(EV0::AbstractMatrix, t::DCDPTmpVars, prim::DynamicDrillModel; vftol=VFTOL, maxit0=45, maxit1=20, kwargs...)
    solve_inf_vfit!(EV0, t, prim; maxit=maxit0, vftol=vftol)

    # try-catch loop in case we have insane parameters that have Pr(no action) = 0, producing a singular IminusTEVp matrix.
    converged, iter, bnds = try
        solve_inf_pfit!(EV0, t, prim; maxit=maxit1, vftol=vftol)
    catch
        @warn "pfit threw error. trying vfit."
        converged, iter, bnds = solve_inf_vfit!(EV0, t, prim; maxit=5000, vftol=vftol)
    end
    return converged, iter, bnds
    # return solve_inf_pfit!(EV0, t, prim; maxit=maxit1, vftol=vftol)
end