# 1. The cost function (returns a scalar)
function cost(params)
  … 
  return cost_value
end

# 2. The in-place gradient function (mutates g)
function grad!(g, params)
  …
  g .= ∇C(params)  # MUTATE the gradient holder. 
  return g          # The package ignores the return value, but returning `g` is good Julia convention.
end

# 3. Wrap them into a CostFunction
core_cost = CostFunction(cost, grad!)

# 4. OPTIONAL: Wrap in a TransformedCost to explore relative/log space, 
# bias parameters, or fix certain parameters. 
# If omitted, MDCProblem will automatically use an empty TransformChain.
chain = TransformChain(LogAbsTransform())
transformed_cost = TransformedCost(core_cost, chain)

# 5. Pass to the MDCProblem
sys = MDCProblem(transformed_cost, p0, dp0, mom)

grad_holder = rand(size(p0))
c0 = cost(p0, grad_holder)
# now grad_holder ≈ zeros(size(grad_holder))

# span defines the length of the curve in the negative and positive directions.
# Negative values ⇒ two curves are evolved in parallel, in the directions ±dp0,
# with lengths span[1] and span[2].

sys = MDCProblem(cost, p0, dp0, mom)
# create an object that stores all the curve information

mdc_curves = MDCSolve(sys; span=MDCSpan(-50.0, 100.0))   # solve the ODE that evolves the MD curve
using Plots
plot(mdc_curves)   # custom plotting of the important features of the MD curve.```


`mdc` is a `MDCSolution`. This type comes with functionality for inspecting / interpolating the curve at different points – all hijacked from DifferentialEquations.jl, of course.



## Features

**Code whatever model you like, however you like**  
- MinimallyDisruptiveCurves.jl does not directly interact with your model. You only need to provide a cost function as described above. Hence there is no constraint on the model simulation step inside the cost function. You could even have your model simulation in another language and call it via `pycall.jl` or `ccall.jl`.

**Avoids the curse of dimensionality**  
- Because the package is not based on sampling parameter space, it does not suffer from the curse of dimensionality in models with many parameters.  
- Generating a minimally disruptive curve usually requires somewhere between a few hundred and a few thousand evaluations of your cost function (and its gradient). That is the major computational cost.  
- If you use reverse‑mode automatic differentiation to compute your cost gradient, the computational complexity of gradient evaluation grows favourably with the number of parameters. That is why machine‑learning models with millions of parameters optimise so fast.

**Multithreaded cost‑function summation**  
- Often you want a model to behave one way for one input, and another way for another input. Each case induces a separate cost function. You can sum cost functions with `sum_losses(arr)`, where `arr` is an array of cost functions. The resulting summed cost function runs each component cost on a separate thread, giving an efficient overall cost evaluation.

**Easily manipulate cost functions**  
- Many MinimallyDisruptiveCurves.jl workflows require regularly fixing / freeing cost‑function parameters, or completely reparameterising cost functions.  
- For example, you might be interested in only allowing a few parameters to change, or you might want to bias the curve to align with a particular parameter.  
- The package provides composable `AbstractTransform` objects that you chain together into a `TransformChain` to define a transformation  

  $$ T : \mathbb{R}^N \to \mathbb{R}^M $$  

  from the optimizer space to the physical space. You just wrap your `CostFunction` in a `TransformedCost` along with your chain:

  ```julia
  chain = TransformChain(LogAbsTransform(), ScaleTransform([1.0, 10.0]))
  C_new = TransformedCost(core_cost, chain)

  to obtain a new differentiable cost function \\( \tilde{C} \\) satisfying  

  $$ \tilde{C}(T(p)) = C(p). $$

- Pre‑loaded transforms include  

  - `LogAbsTransform()`: \\( T = \log \circ |\cdot| \\) for relative‑change‑sensitive parameters.  
  - `ScaleTransform(w)`: \\( T : p_i \mapsto w_i\,p_i \\) for biasing towards or away from particular parameters.  
  - `FixedParamsTransform()`: Fixing / freeing (all except / only) named parameters.  

- The source code for each of these is ≤ 10 lines, so you can easily write your own.


## Caveats

**Your cost function needs to be differentiable.**  
- Your cost function is a mathematical expression of what you want your model to do. A little bit of ingenuity can usually find differentiable cost functions that successfully express this
- Otherwise, switch your brain off and take your cost as the L2 deviation from some ideal model output.
- Or finite difference the gradient with `make_fd_differentiable(cost)`, close your eyes, and hope for the best.


**MinimallyDisruptiveCurves.jl generates 1‑D functional relationships (i.e. *curves*).**  
- There may be higher-dimensional functional relationships (i.e. surfaces or hyper-surfaces) in parameter space, over which model behaviour is approximately preserved. You won't automatically find these using MinimallyDisruptiveCurves.jl. But by evolving multiple curves with different initial conditions / at different points in parameter space, you could get somewhere. 

## Related literature

Not convinced this is the right way to do things? [Raman et al. 2017](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.95.032314) provides a literature review of alternative methods and their relative pros/cons.

## Citing, Acknowledgements

**The algorithm**  (minus a few tweaks)

```bibtex
@article{raman2017delineating,
  title   = {Delineating parameter unidentifiabilities in complex models},
  author  = {Raman, Dhruva V and Anderson, James and Papachristodoulou, Antonis},
  journal = {Physical Review E},
  volume  = {95},
  number  = {3},
  pages   = {032314},
  year    = {2017},
  publisher = {APS}
}